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In  recognition  of  the  negative  repercussions  that  soil  erosion  and  sediment  deposition  can  have  on 
military  preparedness  and  the  environment,  the  Strategic  Environmental  Research  and  Development  Program 
(SERDP)  funded  a  research  effort  between  Fiscal  Years  1993  and  1998  to  develop  geographic  information 
system  (GIS)  tools  and  methods  to  enhance  the  accuracy,  spatial  prediction  and  visual  representation  of 
erosion  and  deposition  on  military  lands.  Because  of  the  variability  in  the  quality  and  availability  of  data  to 
support  erosion/deposition  modeling  efforts,  the  project  sought  to  enhance  both  relatively  simple  empirical 
models  and  more  data-intensive  process-based  models.  In  addition,  the  project  addressed  the  demands  for 
long-term  average  annual  estimates  of  erosion/deposition  and  single  event-based  predictions. 

The  Unit  Stream  Power  Erosion  and  Deposition  (USPED)  model  was  developed  as  a  modification  of 
the  Universal  Soil  Loss  Equation  (USLE).  Whereas  the  USLE  is  a  1 -dimensional  empirical  model,  the 
USPED  model  is  2-dimensional,  considering  upslope  contributing  area,  as  well  as  profile  and  tangential 
curvatures.  The  spatial  patterns  of  erosion,  as  predicted  by  the  USPED  model,  are  much  more  accurate  than 
the  traditional  USLE,  particularly  in  areas  of  complex  topography.  In  addition,  whereas  the  USLE  predicts 
only  erosion,  the  USPED  model  predicts  the  spatial  distribution  of  sediment  deposition.  The  CASC2D  model, 
originally  developed  to  simulate  the  volume  and  distribution  of  runoff  for  single  rainfall  events,  was  modified 
to  include  an  upland  erosion  algorithm.  In  tests  at  the  Goodwin  Creek  watershed  in  Mississippi  and  the 
Henson  Creek  watershed  at  Fort  Hood,  Texas,  the  CASC2D-SED  model  produced  remarkably  accurate  runoff 
hydrographs.  Predictions  of  sediment  discharge  were  less  accurate,  but  fell  within  a  reasonable  range  (-50% 
to  +200%).  The  process-based  SIMulated  Water  Erosion  (SIMWE)  model  was  developed  as  a  multi-dimen¬ 
sional  alternative  to  the  Water  Erosion  Prediction  Project  (WEPP)  model  developed  by  the  US  Agricultural 
Research  Service.  It  borrows  from  the  detachment/transport  theory  of  the  WEPP  model,  but  solves  the 
underlying  continuity  equations  by  Green’s  function  Monte  Carlo  methods  to  provide  the  robustness  neces¬ 
sary  for  spatially  variable  conditions  and  high  resolutions.  It  is  a  landscape  scale  bivariate  model  of  erosion 
and  deposition  by  overland  flow  designed  for  spatially  complex  terrain,  soil  and  cover  conditions. 
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EXECUTI  VE  SUMMARY 


Soil  erosion  and  the  consequent  siltation  of  waterways  have  long  been  major  environmental  concerns 
on  military  installations.  In  recognition  of  the  negative  repercussions  that  soil  erosion  and  deposition  have  on 
military  preparedness  and  the  environment,  the  Strategic  Environmental  Research  and  Development  Program 
(SERDP)  funded  this  project  from  Fiscal  Year  1993  to  Fiscal  Year  1998  to  develop  geographic  information 
system  (GIS)  tools  and  methods  to  enhance  the  accuracy  of  soil  erosion  and  deposition  modeling  on  military 
lands  and  to  facilitate  improved  visual  representation  of  model  results.  The  models  developed  and  enhanced 
during  this  effort  provide  significantly  improved  capability  to  estimate  erosion/deposition  potential  as  an 
input  for  choosing  optimal  land  use  management  and  rehabilitation  programs.  More  accurate  modeling  of 
erosion  and  deposition  assists  land  managers  and  trainers  in  optimizing  the  training  schedules,  delineating 
training  areas,  and  monitoring  changes  over  time.  The  models  also  assist  in  maximizing  availability  of 
military  lands  with  minimal  impact  to  natural  resources,  especially  to  soil  and  vegetation.  The  overall  net 
result  of  this  research  is  improved  ability  to  manage  military  lands  and  reduce  land  maintenance  costs. 

To  achieve  maximum  effectiveness  in  large  areas  of  complex  terrain,  integration  of  erosion  and 
deposition  models  with  GIS  is  essential.  Often,  high  resolution  digital  elevation  models  (DEMs)  required  by 
erosion/deposition  models  must  be  interpolated  from  coarser  resolution  DEMs,  scattered  point  data,  or 
topographic  contours.  For  this  project,  a  method  of  interpolation  by  regularized  spline  with  tension  was 
enhanced  for  the  purpose  of  deriving  high  resolution  DEMs.  The  method  supports  the  combination  of  data 
from  various  sources  with  different  accuracies,  with  the  resulting  surfaces  passing  the  closest  to  the  most 
accurate  data  and  allowing  deviation  from  the  less  accurate  data.  When  compared  with  other  methods  of 
interpolation,  the  resulting  DEM  provides  a  much  more  accurate  representation  of  the  actual  terrain. 

The  Universal  Soil  Loss  Equation  (USLE)  and  its  revision  (RUSLE)  are  the  most  widely  accepted 
erosion  models  in  the  world.  They  are  lumped-parameter  semi-empirical  models  developed  to  predict  long¬ 
term  average  annual  soil  erosion  on  agricultural  fields.  The  models  account  for  topography  by  the  use  of  a 
parameter  that  incorporates  only  slope  length  and  steepness.  Such  a  simplistic  approach  cannot  account  for 
convergence  and  divergence  of  slope  or  for  concavities,  convexities  and  other  irregularities  that  affect  erosion 
and  deposition  processes.  With  SERDP  funding,  the  topographic  parameter  of  these  models  was  replaced  with 
an  analog  that  incorporates  upslope  contributing  area  rather  than  slope  length.  With  this  simple  modification, 
the  USLE  better  reflects  the  impact  of  concentrated  flow  on  soil  erosion.  Further  enhancement  of  the  USLE 
resulted  in  the  Unit  Stream  Power  Erosion  and  Deposition  (USPED)  model  which  is  applicable  to  complex 
slope  geometries.  The  USPED  model  is  most  appropriately  applied  in  situations  where  soil  erosion  is  limited 
only  by  the  ability  of  runoff  to  transport  sediment,  whereas  the  USLE  model  is  best  applied  where  erosion  is 
limited  by  the  ability  of  runoff  to  detach  soil  particles.  The  USPED  model  accurately  predicts  greater  erosion 
on  slope  shoulders  than  on  downslope  positions.  Furthermore,  by  measuring  the  change  in  sediment  transport 
capacity  across  a  GIS  grid  cell,  it  also  predicts  sediment  deposition. 

CASC2D  is  a  2-dimensional,  physically  based  rainfall-runoff  model  that  simulates  spatially  variable 
surface  runoff  for  individual  rainfall  events.  As  a  result  of  SERDP  funding,  the  model  was  applied  to  simulate 
watershed  responses  to  military  training  scenarios.  In  addition,  an  upland  erosion  algorithm  (CASC2D-SED) 
was  added  to  the  model  by  merging  the  physically  based  CASC2D  model  with  the  empirical  USLE  model.  In 
tests  at  the  Goodwin  Creek  watershed  in  Mississippi  and  the  Henson  Creek  watershed  at  Fort  Hood,  Texas, 
the  CASC2D-SED  model  produced  remarkably  accurate  runoff  hydrographs.  Predictions  of  sediment  dis¬ 
charge  were  less  accurate,  but  fell  within  a  reasonable  range  (-50%  to  +200%). 

Over  the  past  decade,  several  process-based  models  have  been  developed  with  the  hope  of  completely 
replacing  the  older  empirical  models.  While  these  models  incorporate  the  impact  of  soil,  cover  and  manage- 
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ment  practices  in  great  detail,  the  description  of  topography  is  overly  simplified.  To  overcome  these  short¬ 
comings,  the  SIMulated  Water  Erosion  (SIMWE)  model  was  developed  with  SERDP  funding.  The  SIMWE 
model  is  based  on  the  solution  of  the  continuity  equation  which  describes  the  flow  of  sediment  over  the 
landscape,  depending  on  a  steady  state  water  flow,  detachment  and  transport  capacities,  and  properties  of  soil 
and  cover.  It  is  a  landscape  scale,  bivariate  model  of  erosion  and  deposition  by  overland  flow  designed  for 
spatially  complex  terrain,  soil  and  cover  conditions.  The  underlying  continuity  equations  are  solved  by 
Green’s  function  Monte  Carlo  methods  to  provide  the  robustness  necessary  for  spatially  variable  conditions 
and  high  resolutions.  Tests  show  that  soil  erosion  and  deposition  estimates  provided  by  the  SIMWE  model 
were  spatially  and  quantitatively  accurate.  In  addition,  investigations  conducted  to  determine  if  the  SIMWE 
model  can  be  used  to  predict  the  consequences  of  common  erosion  and  sediment  control  practices  were 
promising. 

In  conjunction  with  the  development  and  enhancement  of  the  erosion  and  deposition  models,  numer¬ 
ous  GIS  tools  were  developed  and  implemented  to  assist  in  the  visualization  of  the  model  results. 
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1.1  Background 

Soil  erosion  and  the  consequent  siltation  of  waterways  have  long  been  major  environmental  concerns 
on  military  installations.  Accelerated  soil  erosion  results  from  and  ultimately  jeopardizes  military  training  and 
testing.  Among  the  research  and  development  priorities  for  Department  of  Defense  lands,  abatement  of  soil 
erosion  ranks  second  only  to  threatened  and  endangered  species  concerns  (Andrulis  Research  Corporation 
1994).  It  has  been  variously  estimated  that  the  cost  to  restore  damaged  Army  lands  to  tolerable  levels  of  soil 
erosion  could  range  from  $100M  to  $200M  per  year  for  up  to  a  decade.  Annual  maintenance  costs  to  keep  soil 
erosion  at  an  acceptable  level  have  been  estimated  around  $40M  using  existing  technology.  In  an  era  of 
declining  budgets,  the  Department  of  Defense  simply  cannot  afford  such  expenditures.  It  is  paramount  that 
more  cost-effective  measures  be  developed  to  plan  and  implement  erosion  control. 

Numerous  erosion  and  sediment  control  technologies  are  readily  available  to  military  installations. 
These  include  revegetation,  construction  of  earthen  features  such  as  sediment  retention  ponds  and  terraces, 
and  the  use  of  a  wide  variety  of  commercially  available  erosion  control  products.  Unfortunately,  it  is  not 
always  clear  which  techniques  to  use,  where  they  should  be  placed,  what  size  they  should  be,  or  when  the 
optimal  time  occurs  for  implementation.  As  a  result,  projects  are  often  over-  or  under-engineered.  In  addition, 
land  managers  often  address  the  symptoms  of  a  problem  (e.g.,  an  area  of  intensive  erosion)  while  failing  to 
consider  the  ultimate  source  of  the  problem  (e.g.,  the  source  of  excessive  runoff).  Cost-effectiveness  of  land 
reclamation  practices  can  be  maximized  as  we  better  understand  and  model  landscape  processes  that  affect 
soil  erosion.  With  adequate  understanding  and  modeling  capability,  it  is  possible  to  intervene  at  the  appropri¬ 
ate  times  and  places  and  with  the  appropriate  techniques  to  achieve  maximum  benefit  with  the  least  expendi¬ 
ture  of  human  energy  and  financial  resources. 

Until  recently,  most  approaches  to  erosion  and  sediment  modeling  relied  on  lumped-parameter  semi- 
empirical  relationships  developed  for  agricultural  lands.  The  most  widely  accepted  and  used  model  is  the 
Universal  Soil  Loss  Equation  (USLE)  (Wischmeier  and  Smith  1978)  and  its  Revised  form  (RUSLE)  (Renard 
et  al.  1997).  These  equations  estimate  long-term  average  annual  soil  loss  as  the  product  of  factors  represent¬ 
ing  rainfall  erosivity,  inherent  soil  erodibility,  the  length  and  steepness  of  slope,  plant  cover  and  conservation 
support  practices.  While  these  models  are  helpful  in  predicting  average  annual  soil  erosion  over  relatively 
homogenous  parcels  of  land,  they  provide  little  insight  into  the  landscape  processes  governing  soil  erosion, 
and  they  are  incapable  of  accurately  predicting  erosion  on  complex  terrain.  In  addition,  they  cannot  predict 
the  extent  or  spatial  distribution  of  sediment  deposition. 

Over  the  last  decade,  there  has  been  a  move  to  replace  the  empirical  models  with  more  complex 
process-based  models  such  as  the  Water  Erosion  Prediction  Project  (WEPP)  (Flanagan  and  Nearing  1995). 
While  still  developmental  and  not  widely  accepted,  these  models  have  improved  the  understanding  of  ero- 
sional  processes  and  provide  the  basis  for  dramatic  improvements  in  the  accuracy  of  erosion  and  sediment 
modeling. 

A  major  drawback  of  existing  empirical  and  process-based  models  has  been  the  1 -dimensional 
approach  used  to  derive  them.  Landscapes  have  generally  been  assumed  to  be  homogenous,  planar  features. 
Average  erosion  rates  have  been  determined  and  assigned  to  entire  hillslopes  and  watersheds,  thus  providing 
no  information  regarding  the  sources  and  sinks  of  eroded  materials.  Alternatively,  complex  landscapes  have 
been  computationally  divided  into  a  series  of  semi-homogenous  planes,  and  erosion  has  been  calculated  for 
each  plane.  Neither  approach  provides  adequate  spatially  distributed  information  on  erosion  and  deposition  to 
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effectively  optimize  erosion  control  efforts.  Only  through  integration  with  a  geographic  information  system 
(GIS)  is  it  possible  to  model  erosion,  sediment  transport  and  sediment  deposition  in  an  environment  that  can 
account  for  heterogenous  landscapes,  i.e.,  convergence,  divergence,  concavity  and  convexity  of  slopes, 
variable  land  uses  and  conditions,  heterogenous  soil  types,  etc. 

Geographic  information  systems  represent  variability  in  terrain  with  digital  elevation  models 
(DEMs),  where  elevation  is  recorded  for  each  pixel  or  grid  cell  in  a  map.  DEMs  are  generally  derived  from 
scattered  elevation  data  or  topographic  maps.  Interpolation  of  elevation  data  from  known  points  to  all  pixels 
in  a  map  is  generally  required.  While  numerous  methods  are  available,  some  are  more  accurate  than  others. 
Given  the  importance  of  topography  in  erosion  modeling,  the  best  possible  method  of  interpolation  should  be 
used. 


With  an  adequate  DEM,  various  aspects  of  the  hydrologic  cycle  can  be  modeled.  These  include 
rainfall  infiltration,  runoff,  soil  detachment,  sediment  transport  and  sediment  deposition.  Military  land 
managers  are  required  to  make  a  variety  of  decisions  related  to  short-  and  long-term  planning  of  military 
training  and  testing  operations,  and  the  need  for,  timing  and  placement  of  land  reclamation  activities.  The 
most  appropriate  hydrologic  model  for  one  application  may  be  unsuitable  or  overkill  for  another  application. 
Therefore,  a  variety  of  models  should  be  available  to  meet  the  various  data  requirements.  In  addition,  there  is 
a  great  need  to  develop  the  capability  within  these  models  to  predict  the  hydrologic  consequences  of  various 
land  reclamation  practices. 

The  phrase  “a  picture  is  worth  a  thousand  words”  holds  true  for  hydrologic  models.  No  amount  of 
text  or  numerical  output  can  replace  the  value  of  a  picture  that  illustrates  the  spatial  distribution  of  erosion 
and  deposition.  Such  depictions,  as  well  as  3-dimensional  visualization  of  soil  profiles,  watershed  cross- 
sections,  etc.  can  be  invaluable  for  planning  and  placing  erosion  and  sediment  control  projects.  While  multi¬ 
dimensional  visualization  is  common  for  many  applications,  it  has  been  problematic  in  natural  resource  GIS 
applications  where  data  at  a  variety  of  scales  may  need  to  be  combined  (e.g. ,  terrain  surface  features  are 
generally  mapped  at  a  scale  of  meters  while  soil  profile  features  are  measured  in  centimeters). 


1.2  Objectives 

Given  the  deficiencies  identified  in  erosion  and  deposition  models,  as  well  as  the  methodologies  to 

process  digital  elevation  data  and  graphically  represent  the  spatial  and  temporal  distribution  of  natural 

resources  data  to  support  erosion/deposition  models,  the  following  objectives  were  established: 

1 .  Develop  multivariate  spline  interpolation  methods  and  topographic  analysis  tools  to  support  terrain 
modeling  and  processing  of  field  data  (Section  2). 

2.  Enhance  visualization  techniques  supporting  the  design  and  communication  of  dynamic  erosion  and 
sediment  transport  model  results  (Section  3). 

3.  Based  on  upslope  contributing  area  and  the  unit  stream  power  theory,  modify  the  Universal  Soil  Loss 
Equation  (USLE)  to  improve  prediction  of  erosion,  add  prediction  of  deposition,  and  allow  application  of 
the  model  in  complex  topography  (Section  4). 

4.  Add  an  upland  erosion  module  (CASC2D-SED)  to  the  CASC2D  rainfall-runoff  model  to  allow  event- 
based  erosion/deposition  prediction  (Section  5). 

5.  Develop  the  SIMulated  Water  Erosion  (SIMWE)  model  as  a  multi-dimensional  application  of  the  detach¬ 
ment/transport  capacity  theory  approach  to  erosion  and  sediment  prediction  as  contained  in  the  Water 
Erosion  Prediction  Project  (WEPP)  (Section  6). 
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The  U.S.  Army  Corps  of  Engineers  Waterways  Experiment  Station  in  Vicksburg,  Mississippi  was 
responsible  for  Objective  4.  The  remaining  objectives  were  addressed  in  a  collaborative  effort  between  the 
U.S.  Army  Corps  of  Engineers  Construction  Engineering  Laboratory  in  Champaign,  Illinois  and  the  Univer¬ 
sity  of  Illinois. 


1.3  Products 

Each  of  the  objectives  listed  above  was  completed  and  the  corresponding  products  were  placed  in  the 
public  domain.  In  addition  to  this  and  other  periodic  reports  to  the  SERDP  office,  this  project  has  contributed 
significantly  to  the  scientific  literature.  As  a  result  of  efforts  related  to  this  project,  the  authors  have  produced 
over  40  professional  journal  articles,  book  chapters,  conference  proceedings,  professional  presentations  and 
trade  journal  articles.  These  are  listed  in  Appendix  A. 
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2.  DIGITAL  LANDSCAPE 
CHARACTERIZATION  AND  ANALYSIS 


To  model  soil  erosion/deposition  using  a  GIS,  an  adequate  discrete  representation  of  terrain,  soil  and 
cover  properties,  and  climatic  phenomena  must  be  created  in  a  GIS  database.  These  phenomena  are  usually 
measured  at  irregularly  distributed  points  or  digitized  from  isoline  maps.  In  order  to  be  useful  for  most 
erosion/deposition  modeling  applications,  such  data  need  to  be  converted  to  raster  or  grid  format  through  a 
process  of  interpolation.  This  task  can  be  difficult  to  perform  especially  for  geoscientific  data  due  to  a  variety 
of  reasons:  (1)  surfaces  are  complex,  and  the  complexity  is  spatially  heterogeneous,  (2)  the  spatial  distribu¬ 
tion  of  data  points  is  often  heterogeneous  and  inadequate,  (3)  data  sets  are  often  very  large,  thus  creating 
computational  complexity,  and  (4)  data  are  often  noisy. 


2.1  Multivariate  spline  interpolation 

Various  methods  have  been  developed  to  solve  the  problems  of  interpolation  (Mitas  and  Mitasova 
1999),  although  a  general  solution  giving  perfect  results  for  all  data  sets  does  not  exist.  An  illustration  of  the 
results  obtained  by  various  methods  available  in  GIS  is  given  in  Figure  1.  Significant  improvement  in  the 
interpolation  of  digital  elevation  models  and  other  scattered  data  has  been  achieved  by  new  spline  interpola¬ 
tion  and  approximation  methods.  Spline  methods  are  based  on  the  assumption  that  the  approximation  function 
should  pass  as  closely  as  possible  to  the  given  data  points  and  should  be  as  smooth  as  possible.  These  two 
requirements  can  be  combined  into  a  single  condition  of  minimizing  the  deviation  from  the  measured  points 
and  the  smoothness  seminorm  of  the  function.  The  solution  to  this  problem  can  be  expressed  as  a  sum  of  two 
components  (Talmi  and  Gilat  1977):  a  trend  function  and  a  radial  basis  function  with  an  explicit  form  that 
depends  on  the  choice  of  the  smooth  seminorm.  For  our  choice  of  the  smooth  seminorm  (Mitasova  and  Mitas 
1993),  the  trend  function  is  a  constant  and  the  general  form  of  the  radial  basis  function  for  a  t/-variate  case 
(Mitasova  et  al.  1995)  includes  an  incomplete  gamma  function  (Abramowitz  and  Stegun  1964).  We  call  this 
function  Regularized  Spline  with  Tension  (RST).  For  the  special  2D  case,  the  radial  basis  function  can  be 
expressed  through  a  logarithmic  and  exponential  integral  function  (Abramowitz  and  Stegun,  1964).  The 
explicit  form  for  the  function  for  3D  (volume)  and  4D  (volume  in  time)  data  has  been  presented  by  Mitasova 
and  Mitas  (1993)  and  Mitasova  et  al.  (1995).  The  coefficients  of  the  interpolation  function  are  obtained  by 
solving  the  system  of  N  linear  equations,  where  N  is  the  number  of  given  points. 

In  addition  to  high  accuracy  (Mitasova  and  Mitas  1993,  McCauley  and  Engel  1995,  Rohling  et  al.  1998),  the 
RST  has  several  other  useful  properties.  The  generalized  tension  parameter  controls  the  distance  over  which  a 
given  point  influences  the  resulting  surface.  For  the  bivariate  (2D)  case,  tuning  the  tension  can  be  interpreted 
as  tuning  the  character  of  the  resulting  surface  between  a  membrane  and  a  thin  plate.  The  smoothing  param¬ 
eter  can  be  used  to  interpolate  smooth  surfaces  from  noisy  data.  This  parameter  can  be  spatially  variable 
depending  on  the  noise/accuracy  for  each  data  point.  The  proper  choice  of  smoothing  and  tension  parameters 
is  important  for  successful  interpolation  or  approximation.  By  tuning  the  tension  and  smoothing  it  is  possible 
to  minimize  the  overshoots  (Mitasova  and  Mitas  1993),  artificial  pits,  and  the  banding  effect  of  the  elevation 
values  around  contour  lines  as  observed  for  the  less  general  forms  of  splines  such  as  the  thin  plate  spline  or 
thin  plate  spline  with  fixed  tension.  GIS  implementation  of  this  function  provides  several  tools  for  controlling 
the  quality  of  resulting  DEMs,  such  as  computation  of  deviations,  predictive  errors  based  on  ordinary  cross- 
validation  procedures,  and  computation  of  curvatures  useful  for  detecting  artificial  waves  around  contours 
caused  by  improperly  chosen  tension  parameters.  An  optimal  smoothing  parameter  can  be  found  by  ordinary 
or  general  cross-validation  schemes  (Wahba  1990,  Mitasova  et  al.  1995). 


Center  for  Ecological  Management  of  Military  Lands 
Modeling  of  Soil  Erosion 


5 


Figure  1.  Interpolation  of  a  digital  elevation  model  (DEM)  from  scattered  data  points  using  methods  available  in  geo¬ 
graphic  information  systems  (GIS):  (a)  Voronoi  polygons,  (b)  triangulated  irregular  network  (TIN)  based  linear  interpola¬ 
tion,  (c)  inverse  distance  weighting,  (d)  ordinary  kriging,  (e)  spline  with  tension  and  stream  reinforcement  with  Arclnfo 
TOPOGRID,  and  (f)  regularized  spline  with  tension  and  smoothing.  Data  are  from  the  Scheyem  Experimental  Farm, 
Germany  courtesy  of  Dr.  Karl  Auerswald. 
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Due  to  the  solution  of  systems  of  linear  equations,  the  computational  demands  for  the  presented 
method  are  proportional  to  N3  and  their  application  to  large  data  sets  becomes  problematic.  A  solution  to 
interpolation  of  large  data  sets  using  an  explicit  interpolation  function  was  proposed  by  Franke  (1982)  and 
further  developed  and  applied  in  geosciences  by  Mitasova  and  Mitas  (1993).  The  approach  is  based  on  the 
division  of  a  given  region  into  rectangular  segments  and  the  computation  of  an  interpolation  function  for  each 
segment  using  the  data  points  from  this  segment  and  from  its  neighborhood.  Because  the  method  works  with 
equally  sized  segments,  the  approach  is  not  very  efficient  for  strongly  heterogeneous  data  like  digitized 
contours  or  clustered  drill  hole  data.  A  more  effective  approach  to  decomposition  of  a  given  region  with 
heterogeneous  spatial  distribution  of  data  into  segments  with  approximately  the  same  number  of  points  for 
each  segment  has  been  designed  with  quadtrees  (Mitasova  et  al.  1995).  The  interpolation  function  is  then 
computed  for  each  segment  using  the  points  from  this  segment  and  the  points  located  in  the  window  which 
adjusts  its  size  to  the  density  of  points  in  the  neighborhood  of  each  segment.  The  segmentation  is  sensitive  to 
tension;  for  very  low  tension  parameters  used  in  flat  areas,  the  number  of  points  needed  for  smooth  connec¬ 
tion  of  segments  can  be  very  large.  In  certain  extreme  situations  where  a  flat  segment  is  adjacent  to  a  segment 
with  rapidly  changing  terrain,  additional  smoothing  may  be  necessary.  To  improve  the  performance  of  the 
program  for  such  situations  the  segmentation  procedure  was  further  enhanced  by  using  an  automatically 
adjustable  number  of  data  points  dependent  on  the  size  of  each  segment. 


2.2  Topographic  analysis 

Topographic  parameters  needed  for  erosion  modeling  include  parameters  describing  the  local  geom¬ 
etry  of  the  surface,  i.e.,  slope,  aspect,  curvatures,  upslope  contributing  area,  and  slope  length  (Figure  2).  RST 
has  been  specially  designed  to  meet  the  demands  of  topographic  analysis.  It  has  continuous  derivatives  of  all 
orders  allowing  the  use  of  these  derivatives  to  compute  slope,  aspect  and  curvatures  simultaneously  with 
interpolation. 

2.2.1  Local  geometry  parameters.  Regularized  spline  with  tension  (RST)  was  specially  designed  to  support 
direct  computation  of  topographic  parameters  which  are  functions  of  first  and  second  order  derivatives.  Using 
the  partial  derivatives  of  the  RST  function,  slope  angle  is  computed  as  the  arctangent  of  the  magnitude  of  the 
elevation  surface  gradient,  and  aspect  angle  is  computed  as  arctangent  of  the  direction  of  minus  gradient 
(Mitasova  and  Hofierka  1993).  Computation  of  curvatures  is  more  complicated  because,  in  general,  the 
surface  has  different  curvatures  in  different  directions.  A  determination  of  which  curvatures  are  important  is 
dependent  on  the  type  of  processes  under  study.  For  applications  in  geosciences,  curvature  in  the  primary 
gradient  direction  {profile  curvature)  is  important  because  it  reflects  the  change  in  slope  angle  and  thus 
controls  the  change  of  velocity  of  water  flowing  down  along  the  slope  curve.  The  curvature  in  a  direction 
perpendicular  to  the  gradient  reflects  the  change  in  aspect  angle  and  influences  the  divergence/convergence  of 
water  flow.  This  curvature  is  usually  measured  in  the  horizontal  plane  as  the  curvature  of  contours  and  is 
called  plan  curvature  (Zevenbergen  and  Thome  1987,  Moore  et  al.  1991).  However,  for  the  study  of  flow 
divergence/convergence,  it  is  more  appropriate  to  introduce  a  curvature  measured  in  the  normal  plane  in  the 
direction  perpendicular  to  gradient  (Krcho  1991).  This  curvature  is  called  tangential  curvature  because  the 
direction  perpendicular  to  gradient  is,  in  fact,  the  direction  of  a  tangent  to  a  contour  at  a  given  point.  Equa¬ 
tions  for  these  curvatures  can  be  derived  using  the  general  equation  for  curvature  of  a  plane  section  through  a 
point  on  a  surface  as  described  by  Mitasova  and  Hofierka  (1993).  The  positive  and  negative  values  of  profile 
and  tangential  curvature  can  be  combined  to  define  the  basic  geometric  relief  forms  (Krcho  1991,  Dikau 
1990).  Each  form  has  a  different  type  of  flow.  Convex  and  concave  forms  in  the  gradient  direction  have 
accelerated  and  decelerated  flow,  respectively;  convex  and  concave  forms  in  tangential  direction  have  con¬ 
verging  and  diverging  flow,  respectively. 

2.2.2  Flow-related  parameters.  Derivation  of  flow-related  parameters,  such  as  upslope  contributing  area 
and  slope  length  requires  implementation  of  an  effective  flow-tracing  algorithm.  Standard  algorithms  for  flow 
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Figure  2.Topographic  parameters  describing  the  local  geometry  (slope,  aspect  and  curvatures)  and  flow  over  the  terrain 
surface  (flowpath  length,  upslope  contributing  area  and  flowlines). 
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tracing,  which  use  only  a  limited  number  of  flow  directions  (most  often  8)  from  each  grid  cell,  can  lead  to 
various  unrealistic  situations,  such  as  prevailing  flow  in  the  direction  parallel  to  x  or y  axes  or  diagonals 
(Fairfield  and  Leymarie  1991).  Several  new  algorithms  for  flow  tracing  help  to  overcome  deficiencies  of 
standard  algorithms  by  using  the  random  eight  node  approach  (Fairfield  and  Leymarie  1991),  multiple  nearest 
neighbor  nodes  (Freeman  1991),  and  by  using  360  directions  of  flow  with  the  vector-grid  algorithm  (Mitasova 
and  Hofierka  1993,  Mitasova  et  al.  1996). 

2.2.2.1  Flow  path  length.  Flow  path  length  is  used  in  the  standard  form  of  USLE/RUSLE  and  is  appropriate 
for  hillslopes  with  negligible  water  flow  convergence/divergence.  For  the  computation  of  flow  path  length, 
flow  lines  are  generated  uphill  from  each  grid  cell  in  the  gradient  direction  until  they  reach  a  boundary  line, 
singular  point,  barrier,  or  a  grid  cell  with  the  slope  lower  than  the  specified  minimum.  Flow  path  length  for 
each  grid  cell  is  then  computed  from  the  coordinates  of  points  of  their  intersection  with  grid  cells.  Uphill  flow 
lines  converge  on  ridgetops. 

2.2.2.2  Upslope  contributing  area.  Upslope  contributing  area  is  the  entire  area  that  contributes  runoff  into  a 
given  grid  cell.  The  upslope  contributing  area  for  any  given  grid  cell  is  computed  as  the  sum  of  the  area 
comprised  by  all  grid  cells  which  contribute  runoff  into  the  cell  of  concern  (Moore  et  al.  1991).  This  approxi¬ 
mation  is  acceptable  if  the  DEM  is  interpolated  with  adequate  resolution.  Our  experience,  supported  by 
several  recent  studies  (e.g.,  Zhang  and  Montgomery  1994),  is  that  2  m  to  20  m  resolution  is  appropriate  for 
models  using  upslope  contributing  areas  in  regions  with  complex  terrain.  For  the  computation  of  the  number 
of  cells  draining  into  each  grid  cell,  flow  lines  are  constructed  downhill  from  each  grid  cell  in  the  minus 
gradient  direction,  until  they  reach  a  boundary  line,  singular  point,  barrier  (e.g.,  road),  or  a  cell  with  a  slope 
lower  than  the  specified  minimum.  An  improved  algorithm  for  the  construction  of  flow  lines  based  on  vector- 
grid  approach  (Mitasova  and  Hofierka  1993)  is  used.  Flow  lines  constructed  by  this  algorithm  are  represented 
in  vector  format.  The  points  defining  the  flow  line  are  computed  as  the  points  of  intersection  of  a  line  con¬ 
structed  in  the  flow  direction  given  by  aspect  angle  and  the  grid  cell  edges.  Downhill  flow  lines  merge  in 
valleys  and  can  also  be  used  to  define  the  location  of  channels  (Jenson  and  Domingue  1988). 


2.3  GIS  implementation  of  interpolation  and  topographic  analysis 

2.3.1  Geographic  Resources  Analysis  Support  System  (GRASS).  RST  for  bivariate  interpolation  was  first 
released  with  GRASS4.1  as  the  command  s.surf.tps.  The  original  program  has  been  completely  redesigned 
and  enhanced  for  the  GRASS5.0  produced  by  Baylor  University.  The  new  version  of  the  RST  interpolation 
program  was  renamed  to  s.surf.rst  (Appendix  B).  The  most  important  enhancements  include  (1)  output  of 
floating  point  rasters  representing  elevation,  slope,  aspect,  profile  and  tangential  curvatures  as  well  as  first 
and  second  order  partial  derivatives  of  the  modeled  surface,  (2)  output  of  a  site  file  with  deviations  in  given 
points,  (3)  enhanced  segmentation  including  an  option  to  output  the  quadtree  segments,  (4)  spatially  variable 
smoothing  to  support  the  approximation  of  surfaces  from  data  with  heterogeneous  accuracy,  and  (5)  easier  to 
modify  source  code.  If  a  raster  DEM  is  available,  local  geometry  parameters  (slope,  aspect,  curvatures)  can  be 
computed  in  GRASS5.0  with  r.slope.aspect  (Appendix  B).  The  enhancements  include  support  for  floating 
point,  computation  of  curvatures  and  estimation  of  partial  derivatives.  A  vector-grid  flowtracing  algorithm  for 
computation  of  upslope  contributing  area,  hillslope  length  and  flowlines  was  implemented  as  r.flow  (Appen¬ 
dix  B). 

2.3.2  ARC/INFO,  Arc  View.  Spatial  interpolation  can  be  performed  by  TOPOGRID,  which  is  a  different 
version  of  thin  plate  spline  with  tension,  developed  by  Hutchinson  (1989).  The  DEM  derived  by  TOPOGRID 
should  be  used  with  caution  for  erosion/deposition  modeling,  because  it  tends  to  create  waves  along  contours 
leading  to  artificial  spatial  patterns.  Thin  plate  spline  with  tension  and  regularized  spline  developed  by  Mitas 
and  Mitasova  (1988)  are  available  in  ARCGRID  and  ArcView  Spatial  Analyst  as  the  SPLINE  command. 
These  functions  do  not  include  smoothing  capabilities  or  topographic  analysis.  These  must  be  performed 
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separately  after  the  surface  is  interpolated,  using  the  functions  SLOPE,  ASPECT,  CURVATURE  command 
functions.  FLOWACCUMULATION  and  FLOWLENGTH are  derived  using  functions  based  on  the  D8  algo¬ 
rithm  (Jenson  and  Domingue  1988).  This  algorithm  works  well  for  low  resolution  data,  but  at  high  resolutions 
it  tends  to  create  artificial  flow  patterns  biased  towards  grid  cell  diagonals. 
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Implementation  of  landscape  process  simulation  tools  within  a  GIS  has  stimulated  integration  of  GIS 
and  computer  cartography  with  scientific  visualization  (Brown  and  Gerdes  1992,  Brown  and  Astley  1995, 
Brown  et  al.  1995,  Mitasova  et  al.  1994).  Scientific  visualization  is  used  as  both  a  process  of  research  and 
discovery  and  a  method  of  communicating  measured  or  modeled  geographic  phenomena.  As  a  process  of 
discovery,  the  visualization  process  is  cyclical  in  nature,  with  immediate  interactive  visualizations  feeding  a 
refinement  of  the  model.  We  have  found  that  effective  visualization  of  intermediate  results  during  model 
development  is  vital  for  confirming  modeled  phenomena  with  measured  results,  and  encourages  more  rigor¬ 
ous  testing  and  evaluation  of  the  model.  Combinations  of  various  graphical  representations  of  raster,  vector, 
and  point  data  displayed  simultaneously  allow  researchers  to  study  and  query  all  of  their  spatial  data  in  3D 
space.  At  the  same  time,  visual  analysis  of  data  requires  the  capability  to  distort  this  spatial  data  by  changing 
vertical  scale,  separating  surfaces,  performing  simple  transformations  on  point  or  vector  data  for  scenario 
development,  etc. 

Our  visualization  software  tools  provide  methods  and  tools  for  creating  dynamic  cartographic  models 
representing  landscape  phenomena  and  processes.  These  tools  display  multiple  dynamic  surfaces  and 
isosurfaces,  together  with  draped  raster,  vector  and  point  data  in  an  appropriate  projection  of  3D  space  with 
selectable  parameters  of  the  visualization  environment. 

For  interactive  viewer  positioning,  scaling,  zooming,  etc.,  we  use  custom  graphical  user  interface 
(GUI)  widgets  and  a  “fast  display  mode”  where  only  wire  mesh  representations  of  the  data  are  drawn.  When 
rendering  a  scene,  the  user  may  select  various  preset  resolutions  for  better  control  over  rendering  time.  For 
positioning  we  also  chose  to  use  a  paradigm  of  a  moving  viewer  rather  than  a  moving  object  because  we  think 
it  is  more  intuitive  when  modeling  the  reality  of  generally  immobile  geography.  As  a  result,  the  user  feels  as 
though  they  are  flying  around  the  geographic  space,  rather  than  holding  it  in  their  hand  and  rotating  it. 
Horizons  are  kept  level,  resulting  in  less  disorientation  and  fewer  unnatural  viewing  angles.  To  focus  on  a 
particular  object,  the  user  simply  clicks  on  the  object  to  set  a  new  center  of  view. 

Initially,  tools  were  implemented  on  the  Silicon  Graphics  platform  using  a  proprietary  graphics 
library,  IRIS  GL.  While  the  chosen  interface  language,  Tel,  is  publicly  available  for  use  on  multiple  platforms, 
the  underlying  graphics  routines  needed  to  be  converted  to  OpenGL,  a  widely  available  graphic  standard. 
Most  of  the  tools  have  now  been  converted  to  use  the  OpenGL  graphics  API.  Implementation  resulted  in  the 
creation  of  several  specialized  visualization  applications  and  utilities,  a  multi-purpose  viewer  application, 
several  application  programmer  interface  libraries  (APIs)  (see  Appendix  B),  and  enhancements  to  existing 
open  source  software. 

GRASS  GIS,  as  an  open  system  with  a  full  range  of  GIS  capabilities,  has  provided  a  sound  basis  for 
testbed  development  of  visualization  tools.  Each  GRASS  data  type  (raster,  vector,  and  site)  plus  our  own  3D 
grid  format  may  be  used  for  visualization  in  a  single  3D  space.  In  our  implementation,  there  are  four  object 
types  and  various  ways  to  represent  each. 

Surfaces.  A  surface  requires  at  least  one  raster  data  set  to  represent  topography  and  may  use  other  raster  data 
sets  to  represent  attributes  of  color,  transparency,  masking,  and  shininess.  These  data  sets  may  have  been 
derived  from  vector  ( e.g .,  contour)  or  scattered  point  data  using  tools  within  the  GIS.  Users  are  allowed  to  use 
a  constant  value  in  the  place  of  any  raster  data  set  to  produce,  for  example,  a  flat  surface  for  reference  pur¬ 
poses  or  a  constant  color  surface. 
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Vector  sets.  3D  vector  sets  are  not  currently  supported,  so  in  order  to  display  2D  vector  sets  in  3  dimensions, 
they  must  be  associated  with  one  or  more  surfaces.  The  2D  vector  sets  are  then  draped  over  any  of  these 
surfaces. 

Site  Objects.  Point  data  is  represented  by  3D  solid  glyphs.  Attributes  from  the  database  may  be  used  to  define 
the  color,  size,  and  shape  of  the  glyphs.  2D  site  data  must  be  associated  with  one  or  more  surfaces,  and  3D 
site  data  may  be  associated  with  a  surface  ( e.g .,  sampling  sites  measured  as  depth  below  terrain). 

Volumes.  3D  grid  data  may  be  represented  by  isosurfaces  or  cross  sections  at  user-defined  intervals.  Color  of 
the  isosurfaces  may  be  determined  by  threshold  value  or  by  draping  color  representing  a  second  3D  grid  data 
set  over  the  isosurfaces.  Implementation  and  initial  testing  of  a  3D  grid  data  file  format  for  managing  volu¬ 
metric  spatial  data  was  completed.  The  storage  format  and  programmer’s  interface  routines  we  developed 
allow  random  access  to  compressed  floating  point  double  precision  3D  data  with  caching.  It  is  fully  integrated 
within  the  GIS,  using  the  established  database  hierarchy  for  header  and  data  files.  The  resulting  library  was 
used  to  write  several  utility  applications  such  as  r3. in. ascii,  r3.out.ascii,  r3.in.grid3,  r3.mask,  r3.null,  r3.info, 
and  g3.region.  In  addition,  a  program  r3.mkdspf  reads  the  3D  grid  data  and  creates  a  “display”  file  containing 
geometry  for  drawing  isosurfaces  to  represent  the  data  for  visualization  purposes.  Future  work  will  incorpo¬ 
rate  the  3D  grid  file  and  display  file  formats  for  use  within  visualization  tools.  The  specification  and  design 
documents  for  our  3D  grid  API,  and  the  function  prototypes  may  be  found  at: 

http://www2.gis.uiuc.edu:2280/modviz/htdoc/g3d/specification.html 

http://www2.gis.uiuc.edu:2280/modviz/htdoc/g3d/protos.html 

For  visualization  purposes,  point  data  is  often  used  to  derive  an  artificial  surface  or  a  volume  from 
which  isosurfaces  may  be  calculated;  for  modeling  purposes  such  derived  data  is  used  as  input  for  models  that 
work  with  grids  or  volumes.  When  interpolating  terrain  from  point  data,  it  is  easy  to  see  undershoots  or 
overshoots  by  displaying  measured  sites  as  spheres  and  seeing  how  the  interpolated  surface  passes  over, 
under,  or  through  the  spheres.  Parameters  of  the  model  are  then  adjusted  and  the  surface  is  re-interpolated 
until  satisfactory  results  are  obtained.  Similarly,  a  mis-measured  point  is  easy  to  spot  when  visualized  by  3D 
position  and  color.  While  statistical  tools  are  helpful  in  analyzing  data  validity,  we  find  that  many  erroneous 
data  points  that  are  well  within  reasonable  range  and  may  otherwise  be  considered  valid  are  easy  to  detect 
using  spatial-temporal  visualization  tools. 

Qualitative  pre-selection  of  raw  data  points  can  help  eliminate  data  with  low  certainty.  We  developed 
tools  ( d.siter ,  Appendix  B)  that  allow  selection  and  visualization  of  subsets  of  data  directly  within  the  GIS  in 
the  native  multi-attribute  site  format.  From  data  consisting  of  soil  cores  we  were  able  to  quickly  sort  out  cores 
with  a  particular  range  of  particle  size  at  specific  depth  ranges,  or  desired  color  characteristics,  helping  us  to 
verify  soil  horizons. 

Visualization  tools  have  been  improved  to  allow  multiple  attributes/dimensions  of  point  data  to 
determine  the  shape,  size,  color  and  orientation  of  solid  3D  markers  as  well  as  the  3D  position  (SG3d  Up¬ 
grades,  Appendix  B).  Given  a  value  for  an  attribute  of  the  data  at  a  particular  site,  visualization  is  accom¬ 
plished  by  mapping  the  value  to  one  or  more  characteristics  of  the  marker.  In  the  case  of  color,  three  attributes 
(representing  red,  green,  and  blue  components)  may  be  used  to  map  a  single  characteristic,  color.  For  each 
data  attribute,  the  3D  viewer  allows  a  transformation  consisting  of  an  addvalue  and  a  multiplyvalue  which 
converts  the  data  value  to  visualization  units,  preventing  duplication  of  data  in  the  database  solely  for  the 
purpose  of  visualization  and  allowing  the  researcher  to  more  freely  explore  visualization  combinations. 

Variable  shape  is  accomplished  by  having  an  “alternative”  size  field  which  is  used  to  scale  the  marker 
vertically.  Plans  to  allow  scaling  the  marker  differently  in  each  dimension  according  to  attributes  of  the  site 
data  will  further  enhance  the  variability  of  shape.  Using  time  stamps  representing  data  measurement  times,  it 
is  possible  to  visualize  data  collection  sites  and  measured  values  in  compressed  time,  drawing  the  site  mark¬ 
ers  at  intervals  proportionally  representing  the  time  of  measurement,  and  perhaps  using  derived  data  as  a 
background  surface.  Such  methods  may  highlight,  for  example,  a  measurement  interval  in  which  instruments 
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were  improperly  calibrated,  a  period  of  sparse  measurement,  inadvertent  influence  due  to  spatial  measure¬ 
ment  patterns,  or  dramatic  change  that  would  need  to  be  investigated  further.  In  such  cases,  much  of  the 
dataset  may  be  found  to  be  useful  and  valid  and  only  the  bad  data  may  be  discarded,  where  otherwise  the 
conclusion  may  have  been  that  the  entire  dataset  was  unusable. 

Scripting  is  used  to  create  animations  from  series  of  data,  automatically  loading  the  data  sets  and 
rendering  frames  for  the  animation.  Animation  is  useful  for  representation  of  change  in  time,  change  in  a 
modeling  parameter  (Ehlschlaeger  and  Goodchild  1994),  change  in  viewer  position  (fly-bys)  or  change  in 
visible  data.  A  keyframe  technique  is  used  to  generate  animations  when  there  is  no  change  in  data,  e.g.,  to 
create  fly-bys  or  show  a  series  of  isosurfaces  in  volumetric  data.  (Mitasova  et  al.  1994,  Brown  et  al.  1995) 
(Also  see  nviz  Tutorial  in  Appendix  B  ).  Figure  3a  shows  several  frames  from  an  animation  where  a  fence  cut 
is  moved  through  data  to  better  view  underlying  surface  structure.  We  implemented  animation  capabilities  in 
both  2D  (. xganim ,  Figure  3b  and  Appendix  B,  r.out.mpeg.  Appendix  B)  and  3D  (NVIZ,  Appendix  B)  tools. 

Producing  an  animation  of  the  data  involves  redrawing  the  scene  repeatedly  with  slightly  changing 
data  or  viewpoint,  saving  each  resulting  image,  and  then  combining  them  into  a  movie  using  multimedia 
tools.  Our  scripting  mechanism  records  the  actions  of  the  user  while  also  allowing  the  use  of  loops.  An 
integrated  tool  in  the  GUI  provides  methods  for  building  the  script,  including  functions  such  as  Open  Loop, 
Close  Loop,  Open  File  Sequence  Loop,  and  Close  File  Sequence  Loop,  which  help  in  the  creation  of  anima¬ 
tions  by  allowing  a  sequence  of  actions  to  be  applied  to  a  sequence  of  data  files.  Similarly,  a  File  Sequence 
Tool  allows  iteration  of  data  associated  with  the  raster,  vector,  and  site  maps  for  producing  the  most  common 
type  of  data  visualization  animation,  where  data  has  been  precalculated  to  represent  a  timestep  or  change  in 
modeling  parameter  and  there  will  be  no  change  in  viewer  position.  This  tool  also  makes  it  easy  to  add  a 
somewhat  constant  dataset  such  as  roads  for  better  spatial  reference.  For  more  details  on  the  use  of  anima¬ 
tions  for  erosion  modeling  see  Mitas  et  al.  1997. 

Multiple  surfaces  have  proven  useful  to  visualize  boundaries  of  layers  when  examining  soil  core  data 
and  soil  horizon  cross-sections.  Achieving  intuitive  visual  representations  of  these  horizons  presents  a  techni¬ 
cal  challenge  in  terms  of  dimensional  scale,  as  demonstrated  by  Figure  4.  The  vertical  dimension  is  often 
quite  small  relative  to  the  other  two  (e.g.,  soil  depth  of  a  few  meters  vs.  region  dimensions  of  several  kilome¬ 
ters),  so  is  often  exaggerated  when  a  single  surface  is  displayed.  This  exaggeration  is  usually  adequate  to  add 
relief  to  an  otherwise  featureless  surface,  but  in  order  to  separate  close  stratified  layers,  the  required  exag¬ 
geration  grossly  distorts  the  modeled  layers.  If  vertical  translation  of  a  surface  is  used  to  separate  surfaces 
enough  that  they  may  be  viewed  separately,  intersections  between  horizons  and  relative  distances  are  misrep¬ 
resented.  This  is  unacceptable  since  these  may  actually  be  the  features  we  are  interested  in  viewing.  To  study 
differences  between  two  similar  surfaces,  we  use  a  scaled  difference  approach  where  only  the  spatial  distance 
between  surfaces  is  exaggerated,  maintaining  correct  surface  intersections. 

Querying  a  2D  data  set  displayed  as  a  raster  image  can  be  thought  of  as  a  scale  operation  and  a 
translation  operation.  When  a  user  clicks  on  a  pixel,  the  relative  position  in  the  image  is  scaled  by  the  resolu¬ 
tion  of  a  pixel  and  the  north  and  east  offsets  are  added  to  obtain  the  geographical  position.  When  displaying 
surfaces  in  3D  space  with  perspective,  however,  clicking  on  a  pixel  on  the  image  of  the  surface  really  repre¬ 
sents  a  ray  through  3D  space.  The  point  being  queried  is  the  intersection  of  this  ray  with  the  closest  visible 
and  unmasked  part  of  one  of  the  objects  in  the  display.  One  method  for  3D  querying  provided  by  some 
graphics  libraries  is  a  “selection”  method,  where  objects  are  “redrawn”  without  actually  drawing  to  the 
screen,  and  any  objects  drawn  at  the  quety  point  are  returned  as  the  “selected”  objects.  This  method  is  slow 
and  at  best  the  returned  object  is  a  polygon  rather  than  a  specific  point  on  the  polygon. 

When  evaluating  surfaces  made  up  of  relatively  large  polygons  or  irregular  triangular  meshes,  the 
extent  of  a  queried  polygon  may  overlap  or  encompass  several  measured  points,  or  the  polygon  may  be 
textured  by  points  derived  from  a  related  dataset  consisting  of  higher  resolution  data.  Therefore  a  more 
precise  3D  querying  mechanism  became  necessary. 
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Figure  3.  Animation  methods:  (a)  moving  a  cutting  plane  through  a  group  of  surfaces  to  examine  underlying  spatial 
relationships,  and  (b)  xganim  application  allows  up  to  4  series  of  raster  data  to  be  animated  simultaneously  while 
examining  spatial  and  temporal  relationships  in  the  data. 
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No  exaggeration  yields 
little  information  about 
spatial  differences 
between  surfaces. 


Uniform  10X 
exaggeration  overly 
distorts  surfaces. 


No  exaggeration, 
surfaces  separated  - 
surface  intersections  are 
lost. 


10X  exaggeration  of 
difference  between 
surfaces  yields  better 
visualization  of  relative 
differences. 


Figure  4.  Visualization  of  two  similar  surfaces  using  the  scaled  difference  approach  results  in  better  representation  of 
surface  iterations. 
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To  solve  the  general  problem  of  querying  multi-resolution  data  in  3D  space,  we  require  the  user  to 
specify  the  type  of  data  they  are  querying  (surface,  point,  or  vector)  and  then  use  our  knowledge  of  the 
geometry  of  that  data  to  perform  a  geometric  query  in  3D.  Figure  5  shows  that  using  cutting  planes  can  allow 
the  user  to  query  a  specific  location  on  a  surface  that  is  covered  by  another  surface.  This  specialized 
point-on-surface  algorithm  can  be  outlined  as  follows:  (1)  transform  point  on  view  plane  to  a  view  ray,  (2) 
intersect  view  ray  with  convex  polyhedron  defined  by  the  intersection  of  the  parallelepiped  view  region  with 
any  user  defined  cutting  planes,  (3)  if  ray  enters  this  polyhedron,  trace  ray  to  find  any  intersections  with 
visible  and  unmasked  parts  of  any  surfaces,  and  (4)  choose  closest  intersection  to  viewer  (or  return  an  ordered 
list).  Such  point-on-surface  functionality  is  useful  for  3D  data  querying,  setting  center  of  view  and  setting 
center  of  rotation  for  vector  transformations. 


Figure  5.  Querying  multiple  surfaces.  A  view  of  surfaces  as  seen  by  the  user  (above).  Two  cutting  planes,  cl  and  c2,  are 
used  to  see  lower  surfaces  better.  When  the  user  queries  the  image  at  point  P,  the  view  ray  intersects  with  the  visible 
space  enclosed  by  the  convex  polyhedron  V  at  points  E  and  X.  The  resulting  line  is  then  traced  for  intersections  with 
surfaces  and  the  intersection  nearest  the  viewing  position  is  used  to  query  the  database. 
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4.  ENHANCEMENTS  TO  THE  UNIVERSAL 
SOIL  LOSS  EQUATION 


4.1  Enhanced  computation  of  the  effects  of  topography 

The  Universal  Soil  Loss  Equation  is  an  empirical  model  designed  for  the  computation  of  soil  erosion 
in  agricultural  fields.  This  equation  was  developed  for  fields  with  negligible  curvature  and  no  deposition,  and 
where  erosion  is  limited  primarily  by  the  ability  of  rainfall  and  runoff  to  detach  soil  particles  ( i.e .,  detachment 
limited)  as  opposed  to  situations  where  erosion  is  determined  primarily  by  the  capacity  of  runoff  to  transport 
sediment  (i.e.,  transport  limited).  The  equation  has  the  form 

E  =  RKLSCP  (1) 

where  average  annual  soil  erosion  (E)  is  computed  as  the  product  of  factors  representing  the  erosivity  of 
rainfall  and  runoff  (R),  inherent  soil  erodibility  (K),  slope  length  and  steepness  (LS),  plant  cover  (C)  and 
conservation  support  practices  (P)  (Wischmeier  and  Smith  1978,  Renard  et  al.  1991).  Various  modifications 
of  this  equation  are  often  applied  to  the  estimation  of  soil  loss  using  GIS. 

The  topographic  factor  for  USLE  has  been  recently  improved  by  segmenting  irregular  slopes  to 
incorporate  of  the  influence  of  profile  convexity/concavity  and  by  improving  the  empirical  equations  for  the 
computation  of  LS  factor  (Foster  and  Wischmeier  1974,  Desmet  and  Govers  1996,  Renard  et  al.  1991).  To 
incorporate  the  impact  of  flow  convergence,  the  slope  length  factor  (L)  was  replaced  by  upslope  contributing 
area  (Moore  and  Burch  1996,  Desmet  and  Govers  1996).  A  modified  equation  for  computation  of  the  LS 
factor  in  finite  difference  form  for  applications  in  GIS  was  derived  by  Desmet  and  Govers  (1996).  A  simpler, 
continuous  form  of  the  equation  for  computation  of  the  LS  factor  at  any  given  point  on  a  hillslope  was 
developed  by  Mitasova  et  al.  ( 1 996)  as  part  of  the  SERDP  effort.  It  has  the  form 

LS(r)  =  (w+1)  [yf(r)/22.13]m  [sinp(r)/0.0896]”  (2) 

where  r  represents  the  point  location  (x,y),  A  is  the  upslope  contributing  area,  p  is  the  slope  angle,  and  m  and 
n  are  constants.  It  has  been  shown  that  the  values  of  m=0.6  and  n=  1.3  give  results  consistent  with  RUSLE 
LS  factor  for  the  slope  lengths  <100m  and  the  slope  angles  <14  deg  with  negligible  tangential  curvature 
(Moore  and  Wilson  1992). 

The  impact  of  replacing  the  slope  length  by  upslope  area  is  illustrated  in  Figure  6  which  shows  that 
upslope  area  better  reflects  the  impact  of  concentrated  flow  on  increased  erosion.  However,  both  the  USLE 
and  RUSLE  consider  erosion  only  along  the  flow  line  without  the  full  influence  of  flow  convergence/diver¬ 
gence,  and  both  the  standard  and  modified  equations  can  be  properly  applied  only  to  areas  experiencing  net 
erosion  (e.g.,  Desmet  and  Govers  1995).  Because  these  equations  are  incapable  of  predicting  sediment 
deposition,  depositional  areas  must  be  excluded  from  the  analysis.  Therefore,  direct  application  of  USLE/ 
RUSLE  to  complex  terrain  within  GIS  is  rather  restricted. 


4.2  The  Unit  Stream  Power-based  Erosion  and  Deposition  (USPED)  model 

As  part  of  the  SERDP  effort,  the  Unit  Stream  Power-based  Erosion  and  Deposition  (USPED)  model 
was  developed  as  a  further  enhancement  to  the  USLE.  It  predicts  the  spatial  distribution  of  erosion  and 
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Figure  6.  Comparison  of  LS-factors  computed  using  hillslope  length  and  upslope  contributing  area. 


deposition  rates  assuming  steady  state  overland  flow  with  uniform  rainfall  excess  conditions.  The  equation  is 
applicable  to  situations  where  erosion  is  limited  by  the  ability  of  runoff  to  transport  sediment  (i.e.,  detachment 
capacity  exceeds  transport  capacity).  With  this  formulation,  water  flow  q(r)  and  sediment  flow  qs(r)  are 
represented  by  a  bivariate  vector  field  q(r),  qs(r).  Water  flow  can  be  expressed  as  a  function  of  upslope 
contributing  area  (A) 

|q(r)|=^(r)/(r)  (3) 

where  A  is  the  upslope  contributing  area  and  i  is  rainfall  intensity.  The  model  assumes  that  critical  shear  stress 
is  negligible  and  that  the  sediment  flow  rate  is  equal  to  sediment  transport  capacity  T(r).  Sediment  transport 
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capacity  is,  in  turn,  approximated  by  a  power  function  of  slope  (r),  water  flow  q(r)  and  a  transportability 
coefficient  Kfy)  which  is  dependent  on  soil  and  cover  (Julien  and  Simons  1985).  The  equation  is  computed  as 

|q,.(r)l  =  T(r)  =  /f(r)|q(r)|m  sinfXr)”  (4) 

where  m,  n  are  constants  that  depend  on  the  type  of  flow  and  soil  properties.  For  situations  where  rill  erosion 
dominates,  these  parameters  are  usually  set  to  m=  1.6  and  «=  1.3;  where  sheet  erosion  prevails,  they  are  set  to 
m=n=\.0  (Moore  and  Wilson  1992,  Foster  1994).  The  transportability  coefficient  Kt  is  assumed  to  be  approxi¬ 
mately  equivalent  to  the  product  of  the  soil  credibility  (K),  cover  (C)  and  conservation  support  practice  (P) 
parameters  of  the  USLE;  the  rainfall  coefficient  (im)  is  assumed  to  approximate  the  rainfall  (R)  parameter. 
Hence  the  equation  can  be  rewritten  as 


T=  im  K  [^(sinp)"]  CP  (5) 

where  LS  =  Am  (sinP”).  It  is  important  to  note  the  difference  in  the  computation  of  LS  in  this  formulation  as 
compared  to  the  modified  USLE  discussed  earlier  in  this  section. 

Net  erosion/deposition  ( ED )  is  estimated  as  the  divergence  of  the  sediment  flow 

ED( r)  =  div  qs(r)  =  KCP{[grad  h{ r)  •  s(r)]  sin(3(r)  -  h(r)  [fyr)  +  fyr)]  }  (6) 

where  h{ r)  is  water  depth,  .v(r)  is  the  unit  vector  in  the  steepest  slope  direction,  k  (r)  is  the  profile  curvature 
(terrain  curvature  in  the  direction  of  the  steepest  slope),  and  k/y)  is  the  tangential  curvature  (curvature  in  the 
direction  perpendicular  to  the  profile  curvature).  The  topographic  parameters  s(r),  k  ( r),  kfy)  are  computed 
from  the  first  and  second  order  derivatives  of  terrain  surface  approximated  by  the  RST  (Mitasova  and  Mitas 
1993,  Mitasova  and  Hofierka  1993,  Krcho  1991,  see  also  Section  2.1). 

The  spatial  distribution  of  erosion/deposition  is  controlled  by  the  change  in  the  overland  flow  depth 
(first  term)  and  by  the  local  geometry  of  terrain  (second  term),  including  both  profile  and  tangential  curva¬ 
tures.  The  local  acceleration  of  flow  in  both  the  gradient  and  tangential  directions  (related  to  the  profile  and 
tangential  curvatures)  play  equally  important  roles  in  spatial  distribution  of  erosion/deposition.  The  impact  of 
the  tangential  curvature  kt(r)  is  therefore  twofold.  First,  kt( r)  influences  the  water  depth  h{ r)  through  its 
control  of  water  flow  convergence/divergence,  with  tangential  concavity  leading  to  rapid  increase  in  water 
depth  and  an  increase  in  the  potential  for  erosion.  Second,  kjj)  causes  a  local  change  in  sediment  flow 
velocity  which  for  tangential  concavity  has  an  opposite  effect  (reduction  in  sediment  transport),  thus  creating 
potential  for  deposition.  The  interplay  between  the  magnitude  of  water  flow  change  and  both  terrain  curva¬ 
tures  reflected  in  the  bivariate  formulation  therefore  determines  whether  erosion  or  deposition  will  occur.  The 
bivariate  solution  provides  a  sound  theoretical  explanation  for  the  results  of  field  experiments  reported  by 
several  authors  ( e.g .,  Busacca  etal.  1993,  Quine  etal.  1994,  Sutherland  1991)  where  the  highest  erosion  was 
observed  on  divergent  shoulder  elements  and  deposition  on  convergent  footslope  elements. 

A  functionally  equivalent  but  computationally  simpler  formulation  of  the  equation  for  net  erosion/ 
deposition  (ED)  is 


ED  =  div(T  •  s)  = 


d(T cos a) 
dx 


+ 


d(T sin  a) 
dy 


(7) 


where  a  is  the  aspect  or  direction  of  the  steepest  slope,  reported  in  degrees. 

GIS  implementation  of  the  USPED  equation  is  rather  simple,  using  the  tools  developed  for  this 
project  and  implemented  in  GRASS5.0.  A  step  by  step  tutorial  is  available  at  www2.gis.uiuc.edu:2280/ 
modviz/erosion/usped.html.  It  is  important  to  note  that  the  algorithms  available  in  ArcView  and  ArcGrid  for 
interpolation  and  topographic  analysis  are  less  sophisticated  that  those  in  GRASS,  therefore  to  obtain  accept¬ 
able  result  more  attention  should  be  paid  to  the  proper  selection  of  resolution  and  to  the  quality  of  DEM. 
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The  CASC2D  model  was  originally  developed  as  a  2-dimensional  (2D)  physically  based  rainfall- 
runoff  model  to  simulate  spatially  variable  surface  runoff  from  single  rainfall  events  (Julien  et  al.  1995).  It  is 
fully  integrated  with  GIS.  The  model  uses  the  Green  and  Ampt  method  to  model  infiltration,  a  2D  diffusive 
wave  formulation  of  the  Saint  Venant  equations  for  overland  flow,  and  a  ID  solution  of  the  diffusive  wave 
formulation  for  channel  flow  routing.  Outputs  include  runoff  hydrographs  and  maps  of  infiltration  depth, 
surface  moisture,  surface  runoff  depth,  channel  runoff  depth,  and  rate  of  infiltration. 


With  SERDP  funding  the  CASC2D  model  was  modified  to  include  an  upland  erosion  algorithm.  The 
approach  is  based  on  a  sediment  transport  equation  developed  by  Kilinc  and  Richardson  (1973)  who  experi¬ 
mentally  examined  soil  erosion  from  overland  flow  generated  by  simulated  rainfall.  They  developed  a  sedi¬ 
ment  transport  equation  for  sheet  and  rill  erosion  for  bare  sandy  soil: 

qs  =  25500  qo2035  So' 664  (  8) 

where  unit  sediment  discharge  (qs)  is  calculated  in  the  units  of  (tons/m  s)  as  a  function  of  runoff  discharge  (qo) 
and  bed  slope  (So).  This  equation  was  further  modified  by  incorporating  parameters  from  the  USLE  to 
accomodate  various  soil  types,  vegetation  cropping  management  factors,  and  conservation  practices  (Julien  et 
al.  1995): 


<?,=  25500 q2fxS!;m^CP  (9) 

The  numerical  formulation  for  surface  runoff  calculations  stems  from  CASC2D  (Julien  et  al.  1995). 
The  algorithm  for  the  continuity  equation  on  elements  (j,k)  is: 


ht+At(j,k)  =  ht(j,k)  +  ieAt 


t  qtyijJ+D-qtyij-lj) 

~[  pp  +  w 


]A t  (10) 


where  h‘(j,k)  and  ht+l(j,k)  denote  flow  depths  at  the  element  (j,k)  at  times  t  and  t+ 1,  respectively;  ie  is  the 
average  excess  rainfall  rate  over  one  time  step  beginning  from  time  t;  q‘x(k,k+l)  and  q‘x(k-l,k)  describe  unit 
flow  rates  in  the  x-direction  at  time  t,  from  (j,k)  to  (j,k+l),  and  from  0,k-l)  to  (j,k)  consecutively;  likewise  q* 
(j  j+1),  q‘  (j-lj)  denotes  unit  flow  rates  in  the  y-direction  at  time  t,  from  (j,k)  to  (j+l,k),  and  from  (j-l,k)  to 
(j,k)  respectively;  and  W  is  the  grid  size. 

The  momentum  equations  in  the  x  and  y  directions  are  solved  using  the  diffusive  wave  approxima¬ 
tion.  In  the  x-direction,  the  friction  slope  for  the  diffusive  wave  approximation  is  computed  as: 


Stfx(k-U)=Sox(k-lk) 


w 


(11) 
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in  which  the  bed  slope  is  given  by: 


SJk-l,k) 


E(J,k-\)-EU,k) 

W 


(12) 


where  E  represents  the  ground  surface  elevation  of  the  element.  Where  the  friction  slope  is  greater  than  or 
equal  to  zero,  the  calculated  unit  discharge  and  unit  sediment  discharge  for  turbulent  flow  are  given  by: 

q‘x(k-lk)=  )  ~[h'U,k- 1)]3  (13) 

n(j,k- 1) 


and 


q^k-Xk)=e'mlq[(k-Ww%fi-W'm^5CP  (14) 

respectively,  where  qx  implies  unit  discharge  and  qsx  implies  unit  sediment  discharge  in  the  x-direction  at  time 
t  from  (j,k-l)  to  (j,k). 

Where  the  friction  slope  is  less  than  zero,  unit  discharge  and  unit  sediment  discharge  are  calculated 
as: 


q'x(k-lk)=-^-lh'U,k)f 

n(j,k) 

and 


(15) 


=^'in\[(k-UfmiS0,(k-U)'m-~CP  (16) 

respectively,  thus  producing  negative  values  and  implying  that  the  flow  direction  is  actually  from  (j,k)  to  (j,k- 

1) 

The  unit  discharge  and  unit  sediment  discharge  in  the  y-direction  are  similarly  calculated  based  on  the 
sign  of  the  friction  slope  in  the  y-direction.  Once  the  direction  of  flow  and  the  unit  sediment  discharge  have 
been  computed,  the  upland  erosion  is  broken  down  into  three  size  fractions  (sand,  silt,  and  clay)  and  routed 
based  upon  how  much  sediment  is  in  suspension,  previous  deposition,  and  how  much  sediment  has  been 
eroded  from  the  soil  surface.  In  determining  how  much  sediment  is  transported  from  the  outgoing  cell,  the 
model  first  gives  priority  to  the  volume  of  sediment  in  suspension,  secondly  to  the  volume  of  sediment  in 
previous  deposition,  and  lastly  the  remaining  volume  of  sediment,  is  eroded  from  the  soil  surface.  In  order  to 
determine  how  much  sediment  stays  in  suspension  and  how  much  is  deposited  on  the  receiving  cell,  the  trap 
efficiency  for  each  size  fraction  is  computed  using: 


T 


Ei 


-XflJ, 

=  1-  e  hV 


(17) 
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where, 

TEi  =  trap  efficiency  for  each  size  fraction 

X  =  longitudinal  length 

(0i  =  fall  velocity  for  each  size  fraction 

h  =  flow  depth 

V  =  flow  velocity 

The  trap  efficiency  indicates  how  much  sediment  deposit  is  on  the  receiving  cell  for  each  size  fraction,  thus 
the  remaining  volume  of  sediment  (l-TEi)  stays  in  suspension  on  the  receiving  cell. 

The  CASC2D-SED  model  routes  water  and  sediment  from  the  upland  areas  to  the  watershed  outlet. 
Sediment  transport  in  channels  for  single  storm-events  assumes  that  the  change  in  channel  bed  elevation  and 
bank  erosion  processes  are  small  compared  to  upland  erosion  processes.  The  model  keeps  track  of  the  time 
changes  in  the  following  parameters:  rainfall  distribution,  cumulative  infiltration  depth,  surface  runoff  depth, 
suspended  sediment  volume,  sediment  flux,  and  net  aggradation/degradation  for  each  pixel. 
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Over  the  past  decade,  several  process-based  models  have  been  developed  with  the  hope  of  completely 
replacing  the  older  empirical  models.  While  these  models  incorporate  the  impact  of  soil,  cover  and  manage¬ 
ment  practices  in  great  detail,  the  description  of  topography  is  overly  simplified.  To  overcome  these  short¬ 
comings,  the  SIMulated  Water  Erosion  (SIMWE)  model  was  developed  with  SERDP  funding. 


6.1  Model  formulation 

The  methodological  framework  for  the  simulation  of  erosion/deposition  processes  within  the  SIMWE 
model  (SIMulation  of  Water  Erosion)  is  based  on  the  description  of  water  flow  and  sediment  transport 
processes  by  first  principles  equations,  a  concept  outlined  previously,  most  often  for  a  one  dimensional  case 
(Foster  and  Meyer  1972,  Bennet  1974).  Within  our  approach,  inputs  and  outputs  of  the  simulations  are 
represented  by  multivariate  functions  (scalar  or  vector  fields)  rather  than  homogeneous  hillslope  segments. 
The  underlying  continuity  equations  are  solved  by  a  Green’s  function  Monte  Carlo  method  to  provide  the 
robustness  necessary  for  spatially  variable  conditions  and  high  resolutions.  Detailed  descriptions  of  equations 
used  in  this  model  are  given  by  Mitas  and  Mitasova  ( 1 998). 

Two-dimensional  shallow  water  flow  is  described  by  the  bivariate  form  of  the  Saint  Venant  equations 
(Julien  et  al.  1995)  where  the  continuity  of  water  flow  relationship  is  coupled  with  the  momentum  conserva¬ 
tion  equation,  and  the  hydraulic  radius  is  approximated  by  the  normal  flow  depth.  The  system  of  equations  is 
closed  using  the  Manning’s  relationship.  For  this  effort,  we  assume  that  the  solution  of  continuity  and  mo¬ 
mentum  equations  for  a  steady  state  provides  an  adequate  estimate  of  overland  flow  for  land  management 
applications  (Flanagan  and  Nearing  1995).  In  addition,  we  assume  that  the  flow  is  close  to  the  kinematic  wave 
approximation,  but  we  include  a  diffusion-like  term  to  include  the  impact  of  diffusive  wave  effects.  The 
incorporation  of  diffusion  in  water  flow  simulation  is  not  new;  similar  terms  have  been  obtained  in  deriva¬ 
tions  of  diffusion-advection  equations  for  overland  flow  (e.g. .  Lettenmeier  and  Wood  1992).  In  our  reformula¬ 
tion,  we  simplify  the  diffusion  coefficient  to  a  constant  and  we  use  a  modified  diffusion  term.  The  diffusion 
constant  that  we  use  is  rather  small  (approximately  one  order  of  magnitude  smaller  than  the  reciprocal  of 
Manning’s  coefficient)  and  therefore  the  resulting  flow  is  close  to  the  kinematic  regime.  However,  the  diffu¬ 
sion  term  improves  the  kinematic  solution  by  overcoming  small  shallow  pits  common  in  digital  elevation 
models  (DEM)  and  by  smoothing  the  flow  over  slope  discontinuities  or  abrupt  changes  in  Manning’s  coeffi¬ 
cient  (e.g.,  due  to  roads  or  other  anthropogenic  changes  in  elevations  or  cover). 

The  basic  relationship  describing  sediment  transport  by  overland  flow  is  the  continuity  of  sediment 
mass,  which  relates  the  change  in  sediment  storage  over  time,  and  the  change  in  sediment  flow  rate  along  the 
hillslope  to  effective  sources  and  sinks  (Haan  et  al.  1994,  Govindaraju  and  Kawas  1991,  Foster  and  Meyer 
1972,  Bennet  1974).  The  sediment  flow  rate  is  a  function  of  water  flow  and  sediment  concentration.  For 
shallow,  gradually  varied  flow,  the  storage  term  can  be  neglected,  leading  to  a  steady  state  form  of  the  conti¬ 
nuity  equation.  The  sources/sinks  term  is  derived  from  the  assumption  that  the  erosion  and  deposition  rates 
ED(r)  are  proportional  to  the  difference  between  the  sediment  transport  capacity  T(r)  and  the  actual  sediment 
flow  rate  |qs(r)|(Foster  and  Meyer  1 972) 

ED(r)  =  c(r)  [T(r)  -  |qs(r)|]  (18) 
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with  the  first  order  reaction  term  c(r)  dependent  on  soil  and  cover  properties.  The  c(r)  is  obtained  from  the 
relationship  (Foster  and  Meyer  1972) 

ED(r)/Dm(r)  +  |qs(r)|/7(r)  =  1  (19) 

which  states  that  the  ratio  of  erosion/deposition  rate  ED(x)  to  the  detachment  capacity  Dm( r)  plus  the  ratio  of 
the  sediment  flow  to  the  sediment  transport  capacity  is  a  conserved  quantity  (unity).  The  sediment  transport 
capacity  and  detachment  capacity  represent  maximum  potential  sediment  flow  rate  and  maximum  potential 
detachment  rate,  respectively,  and  are  functions  of  a  shear  stress  /( r)  (Foster  and  Meyer  1972) 

t(r)  =  w  h  (r)  sinp(r)  (20) 

where  w  is  the  hydrostatic  pressure  of  water  with  the  unit  height  (h).  Then  the  simplified  equation  for  trans¬ 
port  capacity  is 


T(r)  =  Kl(r)t(ry  (21) 

where  Kt( r)  is  the  effective  soil  transportability  coefficient  and  p  is  exponent  dependent  on  the  type  of  flow. 
Detachment  capacity  estimated  as  a  function  of  shear  stress  can  be  expressed  as 

Dm(r)  =  KJj)  [t(ryts(r)y  (22) 

where  K/r)  is  the  effective  erodibility  (detachment  capacity  coefficient),  q  is  an  exponent  parameter  and  /  (r) 
is  the  critical  shear  stress.  The  parameters  and  adjustment  factors  for  the  estimation  of  detachment  and 
transport  capacity  are  functions  of  soil  and  cover  properties.  Typical  values  are  being  developed  for  the 
WEPP  model  for  a  wide  range  of  soils,  cover,  agricultural  and  erosion  prevention  practices  (Flanagan  and 
Nearing  1995). 

To  solve  the  bivariate  continuity  equation,  we  introduce  a  small  diffusion  term  (Mitas  and  Mitasova 
1998)  which  represents  local  dispersion  processes  of  the  suspended  flow  (e.g.,  the  impact  of  small,  local  slope 
changes  below  the  DEM  resolution)(Bennet  1974).  On  the  left  hand  side  of  the  equation,  the  first  term 
describes  local  diffusion,  the  second  term  is  a  drift  driven  by  the  water  flow  while  the  third  term  represents  a 
velocity  dependent  ‘potential’ .  The  size  of  the  diffusion  constant  is  about  one  order  of  magnitude  smaller 
than  the  reciprocal  Manning’s  constant  so  that  the  impact  of  the  diffusion  term  was  relatively  small. 

The  equations  can  be  solved  by  projection  methods  (Rouhi  and  Wright  1995).  An  alternative  is  to 
interpret  equations  as  a  representations  of  stochastic  processes  with  diffusion  and  drift  components  (Fokker- 
Planck  equations)  and  then  carry  out  the  actual  simulation  of  the  underlying  process  utilizing  stochastic 
methods  (Gardiner  1985).  This  is  very  similar  to  Monte  Carlo  methods  in  computational  fluid  dynamics  or  to 
quantum  Monte  Carlo  approaches  for  solving  the  Schrodinger  equation  (Hammond  et  al.  1994,  Mitas  1996). 

The  Monte  Carlo  technique  has  several  unique  advantages  which  are  becoming  even  more  important 
due  to  new  developments  in  computer  technology.  Perhaps  one  of  the  most  significant  Monte  Carlo  properties 
is  robustness  which  enables  us  to  solve  the  equations  for  complex  cases,  such  as  discontinuities  in  the  coeffi¬ 
cients  of  differential  operators  (in  our  case,  abrupt  slope  or  cover  changes,  etc).  Also,  rough  solutions  can  be 
estimated  rather  quickly,  allowing  us  to  carry  out  preliminary  quantitative  studies  or  to  rapidly  extract  qualita¬ 
tive  trends  by  parameter  scans.  In  addition,  the  stochastic  methods  are  tailored  to  the  new  generation  of 
computers  (i.e.,  they  provide  scalability  from  a  single  workstation  to  large  parallel  machines  due  to  the 
independence  of  sampling  points).  Therefore,  the  methods  are  useful  both  for  everyday  exploratory  work 
using  a  desktop  computer  and  for  large  cutting-edge  applications  using  high  performance  computing. 


26 


Center  for  Ecological  Management  of  Military  Lands 
Modeling  of  Soil  Erosion 


6.2  Model  properties 

Below  we  analyze  the  role  of  the  SIMWE  model  parameters  for  complex  terrain  assuming  spatially 
uniform  rainfall  excess,  soil  and  cover  properties. 

The  first-order  reaction  coefficient  c(r)  is  related  to  the  soil  detachability  and  transportability  and 
controls  the  spatial  extent  of  deposition.  There  are  two  limiting  cases  of  erosion  and  sediment  transport,  ie., 
detachment  limited  and  sediment  transport  capacity  limited  (Foster  and  Meyer  1972,  Hairsine  and  Rose 
1992).  The  first  case  occurs  when  c(r)  approaches  0,  resulting  in  the  net  erosion  being  approximately  equal  to 
the  detachment  capacity.  Therefore,  for  conditions  when  c(r)«l,  detachment  limited  erosion  prevails.  In  this 
situation,  almost  all  detached  sediment  is  transported  to  the  stream  and  deposition  is  restricted  to  small 
concave  areas  and  channels.  This  case  for  a  ID  hillslope  is  modeled  by  the  Universal  Soil  Loss  Equation 
(USLE).  For  the  second  limiting  case  represented  by  c(r)  >  1,  sediment  flow  approximates  sediment  transport 
capacity.  In  this  case,  net  erosion/deposition  is  modeled  as  a  divergence  of  the  sediment  transport  flow,  and 
the  model  predicts  large  areas  with  deposition.  Such  a  behavior  is  close  to  the  observed  distribution  of 
colluvial  deposits,  suggesting  the  prevailing  influence  of  the  transport  capacity  limited  case  on  long  term 
patterns  of  deposition.  This  case  is  modeled  by  the  USPED  model  (Mitasova  et  al.  1996).  Simulations  in 
which  the  c(r)  values  increase  from  0.01  (fine  soils)  to  10.0  (sandy  soils),  assuming  a  rough  cover  (grass)  with 
n=0.1,  illustrate  that  c(r)  has  a  dramatic  effect  on  the  spatial  distribution  of  erosion  and  deposition  over  the 
landscape  (Mitas  and  Mitasova  1998).  As  c(r)  increases,  depositional  areas  expand  while  the  sediment  flow 
rate  in  the  streams  decreases  until  the  pattern  typical  for  a  transport  capacity  limiting  case  is  reached.  Because 
the  c(r)  is  dependent  on  the  ratio  between  detachment  capacity  and  sediment  transport  capacity,  both  credibil¬ 
ity  and  transportability  coefficients  influence  the  spatial  extent  of  deposition.  However,  as  we  demonstrate  in 
the  next  section,  each  of  these  coefficients  has  a  different  impact  on  the  sediment  loads  in  streams. 

Erodibility,  represented  by  the  detachment  capacity  coefficient  Kf  r),  is  a  measure  of  the  susceptibility 
of  soil  to  detachment  by  water  flow  (Flanagan  and  Nearing  1995).  It  is  often  defined  as  the  increase  in  soil 
detachment  per  unit  increase  in  shear  stress  of  clear  water  flow.  A  change  in  the  erodibility  factor,  while 
keeping  the  other  parameters  constant,  changes  the  ratio  between  transport  capacity  and  detachment  capacity. 
This  leads  to  the  change  in  the  character  of  erosion  process.  Detachment  capacity  limited  erosion  occurs  when 
erodibility  is  significantly  lower  than  transportability;  transport  capacity  limited  erosion  occurs  when  erod¬ 
ibility  is  greater  than  transportability.  Therefore,  a  change  in  erodibility  will  change  the  spatial  pattern  of 
erosion  and  deposition  (Figure  7a-c).  However,  changes  in  erodibility  have  minimal  impact  on  the  amount  of 
sediment  load  in  the  stream.  This  illustrates  the  problems  associated  with  the  use  of  in-stream  sediment  loads 
for  making  land  management  decisions. 

Transportability,  represented  by  the  transport  capacity  coefficient  Kfy),  is  a  measure  of  the  likelihood 
that  soil  particles  will  be  transported  by  water  flow.  It  depends  on  soil  properties,  but  can  be  also  influenced 
by  vegetation.  This  coefficient  is  not  directly  measured  and  provided  in  the  WEPP;  rather  it  is  estimated 
indirectly,  making  the  proper  determination  of  this  parameter  problematic.  However,  the  parameter  can  be 
derived  at  least  for  some  types  of  soils  using  the  published  values  of  first  order  reaction  coefficient  or  using 
the  procedure  suggested  by  Finkner  et  al.  (1989).  Our  simulations  show  that  this  parameter  has  a  profound 
impact  on  the  erosion  process  as  it  influences  both  spatial  distribution  and  magnitude  of  sediment  flow  and 
erosion/deposition  rates.  Recently,  the  importance  of  transport  capacity  for  overland  flow  erosion  processes 
has  been  fully  recognized  and  more  experimental  and  theoretical  work  is  being  performed  (Guy  et  al.  1991, 
Govers  1991,  Nearing  etal.  1997).  Figures  8(a-c)  ^illustrates  how  the  change  in  transportability  changes  the 
erosion  regime  from  detachment  capacity  limited  to  transport  capacity  limited  while  also  reducing  the  magni¬ 
tude  of  erosion  rates. 

Surface  roughness  is  represented  by  the  Manning’s  n( r).  It  has  a  significant  impact  on  the  location 
and  magnitude  of  deposition.  For  the  same  value  of  c  =  1.0,  the  extent  of  deposition  for  smooth  surfaces  (bare 
soil  with  n=0.01)  is  smaller  than  for  rough  surfaces  (grass  with  n=0.1).  This  also  assumes  an  increase  in 
detachability  and  transportability  for  bare  soil  as  compared  to  grass. 
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Figure  7.  Changes  in  sediment  flux  (left)  and  erosion/deposition  (right)  as  a  result  of  a  change  in  erodibility 
(Kd)  when  transportability  (K)  is  held  constant  at  0.3  [(a)  K  =0.003,  c=0.01;  (b)  ^=0.03,  c=0.1;  (c)  Kj=03, 
c=1.0]. 
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Figure  8.  Changes  in  sediment  flux  (left)  and  erosion/deposition  (right)  as  a  result  of  a  change  in  transport¬ 
ability  (K)  when  erodibility  (Kj  is  held  constant  at  0.003  [(a)  AT =0.3,  c=0.01;  (b)  .£=0.03,  c=0.1;  (c) 
£=0.003,  c=1.0)]. _ 
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Rainfall  excess  /( r)  is  estimated  as  rainfall  intensity  minus  infiltration  rate.  Rainfall  excess  can  be 
estimated  per  unit  time,  as  in  a  dynamic  situation  where  rainfall  intensity  and  infiltration  rate  vary  with  time, 
or  as  a  constant  where  rainfall  intensity  and  infiltration  are  assumed  to  be  in  a  steady  state.  Rainfall  excess 
influences  the  magnitude  of  erosion/deposition  rates.  When  rainfall  excess  increases,  erosion  and  deposition 
rates  also  increase.  However  the  spatial  pattern  of  erosion  and  deposition  does  not  change. 

Critical  shear  stress  t(r)  represents  the  soil’s  resistance  to  the  shearing  forces  (shear  stress)  of 
flowing  water.  It  depends  on  soil  and  cover  properties.  Typical  values  for  many  soils  are  available  (Flanagan 
and  Nearing  1995).  If  the  shear  stress  of  flowing  water  at  the  given  location  is  lower  than  the  critical  shear 
stress  of  the  soil,  no  soil  is  detached.  This  parameter,  therefore,  has  an  impact  on  erosion/deposition  patterns. 
High  values  reduce  the  spatial  extent  of  erosion.  On  the  other  hand,  it  can  also  increase  the  magnitude  of 
erosion  rates  on  steeper  hillslopes  and  in  areas  with  concentrated  flow  due  to  the  fact  that  clean  water  has 
higher  potential  to  transport  the  sediment. 

It  is  important  to  note  that  the  parameters  do  not  act  independently.  They  are  interrelated,  and  it  is 
their  interaction  that  controls  the  pattern  and  magnitude  of  erosion.  For  example,  the  growth  of  vegetation 
reduces  Kd  and  Kt,  and  increases  n.  The  resulting  erosion/deposition  pattern  depends  on  the  interaction 
between  the  rates  of  this  change.  If  both  Kd  and  Kt  change  at  the  same  rate,  the  spatial  distribution  of  erosion/ 
deposition  stays  the  same  and  only  the  magnitude  of  rates  changes.  If  vegetation  growth  reduces  Kd  faster 
than  Kt,  the  erosion/deposition  pattern  will  change  from  transport  capacity  limited  towards  detachment 
limited.  The  relationships  between  these  parameters  is  an  area  ripe  for  additional  research.  There  is  a  lack  of 
systematic  experimental  and  theoretical  modeling  work  in  this  area. 

The  present  analysis  shows  that  for  uniform  soil  and  cover  there  is  a  generalized  pattern  of  erosion 
and  deposition.  High  erosion  risk  areas  are  located  on  upper  convex  parts  of  hillslopes  and  in  hollows  and 
centers  of  valleys  with  concentrated  flow.  Deposition  occurs  in  valleys  and  concave  lower  parts  of  hillslopes. 
By  varying  soil  and  cover  conditions  uniformly  across  a  hillslope,  the  position  where  deposition  starts  can 
move  up  or  down  slope.  Spatial  variability  in  soil  and  cover  has  a  great  impact  on  the  basic  pattern  of  erosion/ 
deposition.  Depending  on  the  spatial  distribution  of  soil  and  cover  parameters,  the  distribution  of  erosion  and 
deposition  may  change.  In  addition,  the  overall  soil  loss  and  sediment  loads  in  streams  may  be  altered. 
Transition  zones  between  land  covers  types  (e.g.,  bare  soil  and  dense  grass)  may  cause  abrupt  changes  in  flow 
velocities,  as  well  as  transport  and  detachment  capacities,  creating  effects  important  for  erosion  prevention. 

Figure  9a  illustrates  this  phenomenon  at  the  Scheyem  Experimental  Farm  in  Germany.  We  performed 
a  simulation  with  the  parameters  set  for  dense  grass  in  the  meadow  area  (n=0.1,  K  =  Kd  =  0.0003 )  and  bare 
soil  in  the  arable  area  ( n-0.01 ,  Kt  =  Kd  =  0.03),  and  compared  the  results  with  the  observed  spatial  patterns  of 
erosion/deposition.  The  highest  rates  of  both  the  net  erosion  and  the  net  deposition  are  predicted  in  valleys 
with  high  concentrated  sediment  flow.  Field  measurements  confirm  that  this  area  has  the  thickest  layers  of 
colluvial  deposits  with  large  linear  erosion  features  observed  after  an  unusually  strong  storm  (Auerswald  et 
al.  1996).  Relatively  high  erosion  rates  were  also  predicted  on  bare  upper  parts  of  hillslopes.  These  predicted 
soil  losses  correlate  well  with  the  loss  of  radio-tracers  and  declines  in  crop  productivity  (Figure  10).  High 
erosion  was  also  predicted  for  narrow  strips  below  the  grass  areas,  where  water  accelerates  after  depositing 
the  sediment  and  leaving  the  grass.  This  leads  to  an  increase  in  the  difference  between  the  sediment  transport 
capacity  and  the  actual  sediment  flow,  creating  the  conditions  for  higher  net  erosion. 

Deposition  is  predicted  at  the  lower,  concave  parts  of  hillslopes  which  is  in  agreement  with  the  spatial 
distributions  of  colluvial  deposits  and  radio-tracers  (Figures  9a  and  10).  The  sediment  flow  rate  decreases 
sharply  and  eroded  materials  are  deposited  at  the  upper  edges  of  grassed  meadows,  thus  explaining  why 
rilling  ends  abruptly  at  the  border  between  the  grassed  and  bare  areas  (Figure  9a).  The  influence  of  grass 
cover  on  flow  velocity,  transport  and  detachment  capacities  also  explains  the  observed  location  of  deposited 
material  on  the  convex,  upper  part  of  the  grassed  hillslopes. 
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Figure  9.  Simulation  of  runoff  depth,  sediment  flow  and  net  erosion/deposition  for  3  different  land  manage¬ 
ment  scenarios  at  the  Scheyem  Experimental  Farm:  (a)  traditional  land  use,  (b)  “best  management  practices,” 
and  (c)  computer  designed  scenario  to  minimize  erosion. 


FigurelO.  Spatial  distribution  of  net  erosion  and  deposition  predicted  by  the  SIMWE  model  and  compared 
with  (a)  rills  observed  after  an  intense  rainstorm  and  (b)  erosion/deposition  estimates  from  radiotracer  data, 
average  crop  yield,  and  measured  colluvial  deposits. 
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7.  APPLICATIONS  AND  CASE  STUDIES 


One  of  the  goals  of  this  project  was  to  develop  methods  and  tools  for  erosion  risk  assessment  and 
erosion  prevention  in  support  of  military  installations.  This  task  poses  a  special  challenge,  because  military 
installations  often  occupy  large  areas  with  terrains  much  more  complex  than  typical  agricultural  regions  for 
which  most  of  the  traditional  erosion  modeling  tools  were  developed.  In  addition,  the  manner  in  which 
military  installations  are  used  often  creates  a  mosaic  of  relatively  well-preserved  natural  areas  intermingled 
with  landscapes  exposed  to  high  intensity  disturbance.  The  principles  of  process-based  erosion  modeling 
developed  for  agriculture  lands  have  to  be  significantly  enhanced  and  new  approaches  have  to  be  developed 
to  meet  this  challenge.  This  section  demonstrates  that  application  of  the  new  technologies  to  real  world 
problems  at  several  sites. 


7.1  The  effects  of  DEM  resolution  (Yakima  Training  Center,  WA) 

In  the  early  stages  of  the  project,  we  tested  the  reliability  of  using  a  standard  U.S.  Geological  Survey 
(USGS)  30m  DEM  to  predict  topographic  potential  for  erosion/deposition  at  the  Yakima  Training  Center, 
Washington  (Mitasova  et  al.  1996).  We  found  that  the  conversion  of  USGS  contour  maps  to  digital  elevation 
models  with  horizontal  and  vertical  resolutions  of  30  m  and  1  m,  respectively,  resulted  in  systematic  banding 
or  artificial  waves  along  the  contours  (Figure  1  la).  To  improve  the  results,  we  used  the  regularized  spline 
with  tension  (RST)  method  (see  Section  2.1)  to  reinterpolate  the  elevations  to  the  new  DEM  with  10m 
horizontal  and  0.01m  vertical  resolution.  Detailed  analysis  of  the  original  and  reinterpolated  DEM,  using  3 
dimensional  shaded  views,  curvatures  and  histograms,  revealed  that  systematic  bands  or  waves  along  the 
contours  were  significantly  reduced  but  not  eliminated  (Mitasova  et  al.  1996).  Profile  curvature  and  histo¬ 
gram  analysis  proved  that  it  was  not  possible  to  remove  the  waves  completely.  While  resampling  and  smooth¬ 
ing  improved  the  results,  the  USPED  model  still  predicted  waves  of  deposition  along  20m  contours  (Figure 
1  lb).  To  confirm  that  the  wave  patterns  were  artificial,  we  extracted  10m  contours  from  the  original  grid 
DEM  and  interpolated  a  new  DEM  using  the  RST  function  with  smoothing  and  tension.  The  new  DEM  had  a 
root  mean  square  error  of  2.3m  which  is  well  within  the  vertical  accuracy  of  the  original  data  (given  as  7m), 
and  the  banding  effect  around  the  contours  disappeared  (Figure  1  lc).  A  comparison  of  the  results  from  the 
original,  reinterpolated  and  recomputed  DEMs  demonstrates  inadequacy  of  traditional  30m  resolution  USGS 
DEMs  for  erosion/deposition  modeling.  In  the  absence  of  higher  resolution  data,  it  is  possible  only  to  com¬ 
pute  the  topographic  potential  for  detachment  limited  erosion  which  does  not  require  the  computation  of 
deposition  and  is  thus  less  sensitive  to  artifacts  in  DEM. 


7.2  Analysis  of  topographic  potential  for  erosion/deposition  at  Ft.  McCoy,  WI 

Topographic  potential  for  soil  detachment,  net  erosion  and  deposition  at  Fort  McCoy  was  estimated 
by  the  enhanced  USLE  and  USPED  models  (Mitasova  et  al.  1996).  Even  with  rather  crude  elevation  data  (30 
m),  it  was  possible  to  identify  some  areas  with  high  topographic  risk  for  soil  erosion  and  deposition  (Figure 
12).  The  plateaus  in  low  elevation  areas  due  to  the  insufficient  vertical  resolution  of  lm  create  an  artificial 
pattern  of  steeper  slopes  along  lm  contours  and  upslope  contributing  areas  computed  from  the  original  30m 
DEM  by  r.flow  are  underestimated.  Pits  and  plateaus  in  the  DEM  cause  problems  in  flowtracing  with  almost 
no  flow  generated  in  lower  elevation  areas  and  disrupted  flow  in  valleys.  The  modified  LS  factor  predicts 
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Figurell.  Net  erosion  and  deposition  predicted  by  the  USPED  model  using  the  original  30  m  DEM  (a),  a  resampled  and 
smoothed  DEM  (b),  and  a  DEM  reinterpolated  from  contours  (c). 


well  some  areas  of  high  potential  for  soil  detachment,  especially  in  hilly  parts  of  the  region.  However,  it  does 
not  predict  the  location  of  areas  with  concentrated  flow  with  high  potential  for  gully  formation.  Also,  in  lower 
elevation  areas  the  artificial  pattern  of  increased  erosion  is  predicted  along  lm  contours.  The  net  erosion/ 
deposition  map  computed  by  the  USPED  model  predicts  high  net  erosion  in  hilly  regions  and  on  slopes  along 
the  streams.  It  also  shows  that  a  significant  portion  of  the  material  eroded  from  hillslopes  is  deposited  before 
it  could  reach  the  main  streams.  However,  this  map  lacks  the  prediction  of  high  erosion  due  to  concentrated 
flow  in  valleys  and  shows  artificial  waves  of  erosion  and  deposition  along  the  lm  contours  in  flatter  areas. 

To  reduce  the  negative  impact  of  the  low  resolution  of  the  available  elevation  data,  the  DEM  was 
reinterpolated  to  10m  horizontal  and  0.01m  vertical  resolution  using  the  RST  method.  Simultaneously  with 
interpolation,  maps  of  slope,  aspect,  profile  and  tangential  curvatures  were  generated.  Steady  state  water  flow 
computed  by  r.flow  from  the  smoothed  10m  DEM  shows  potential  for  channel  formation  in  valleys  and 
predicts  water  flow  also  in  low  elevation  region  (Figure  13).  Flow  in  the  two  main  streams  is  not  adequately 
described  because  of  lack  of  data,  but  it  can  be  incorporated  if  stream  data  are  available.  Modified  LS  factors 
were  computed  from  the  10m  DEM  using  different  exponents  for  the  water  flow  term  ( p=0.6  and  p=l-5). 
Value  of  this  exponent  is  still  a  subject  of  research  and  discussion  in  erosion  research  community.  The  first 
result  for  p=0.6  puts  more  weight  on  the  influence  of  slope  and  theoretically  represents  the  detachment 
limited  erosion  (Figure  13).  The  second  result  with  p=1.5  puts  more  weight  on  the  influence  of  water  flow 
and  theoretically  represents  the  sediment  transport  capacity  limited  erosion  (Figure  13).  Because  the  exponent 
depends  on  the  conditions  of  water  flow  in  a  particular  area  it  should  be  calibrated  to  reflect  the  type  of  flow 
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Figurel2.  Topographic  potential  for  erosion  as  predicted  by  the  enhanced  USLE  (a)  and  erosion/deposition  as  predicted 
by  the  USPED  model  (b)  from  a  30  m  DEM.  The  arrows  in  b  show  artificial  waves  of  deposition  due  to  waves  along  the 
contours  of  the  DEM. 
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Figure  13.  Results  of  erosion  modeling  based  on  a  smoothed  DEM  that  was  resampled  to  10  m  resolution,  (a)  Enhanced 
USLE  LS-factor  with  p=0.6,  representing  topographic  potential  for  detachment  limited  erosion,  (b)  Sediment  transport 
capacity  estimated  as  a  function  of  upslope  area  and  slope  with  p=l  .5.  (c)  Topographic  potential  for  net  erosion  and 
deposition  estimated  as  divergence  of  sediment  flow  with  the  USPED  model  where  the  magnitude  of  sediment  flow  is 
approximated  by  sediment  transport  capacity. 

typical  for  the  modeled  area  and  time  during  the  year.  Topographic  potential  for  net  erosion/deposition  was 
computed  from  the  10m  DEM  using  the  USPED  model.  Similarly  as  for  the  result  of  30m  DEM,  the  model 
shows  high  erosion  in  hilly  area  and  along  main  streams  and  deposition  in  concave  areas.  It  also  indicates 
location  of  areas  with  concentrated  flow  which  can  reach  streams. 


7.3  Impact  of  proposed  land  use  change  at  Camp  Shelby,  MS 

We  used  the  USPED  model  to  evaluate  potential  changes  in  the  spatial  distribution  of  erosion  and 
deposition  at  Camp  Shelby,  Mississippi  resulting  from  anticipated  changes  in  land  use  patterns.  Camp  Shelby 
is  largely  wooded.  In  an  attempt  to  open  more  area  to  military  training,  it  was  proposed  that  large  areas  be 
deforested  and  additional  areas  be  selectively  thinned  to  allow  passage  of  tracked  vehicles.  A  high  resolution 
(5m)  DEM  with  slope  and  aspect  was  created  from  digitized  contours  using  regularized  spline  with  tension 
(see  Section  2.1).  The  upslope  contributing  area  for  each  grid  cell  was  computed  by  the  vector-grid  based 
flow  tracing  algorithm.  A  net  erosion  and  deposition  index  was  computed  using  the  USPED  model  for  the 
current  undisturbed  landuse  and  for  several  land  use  alternatives  with  a  significant  portion  of  the  area  exposed 
to  intensive  landuse  during  training.  One  alternative  is  shown  in  Figure  14.  Calculations  indicate  that  the 
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proposed  scenario,  if  implemented  without  erosion  protection  measures,  will  lead  to  a  3  3 -fold  increase  in 
average  soil  loss  per  acre.  Use  of  protective  stream  buffers  and  exclusion  slopes  steeper  than  10%  from 
training  could  reduce  potential  erosion  by  about  40%,  but  it  would  still  be  22  times  greater  than  before 
clearing  (Figure  14). 


7.4  The  effects  of  scale  on  erosion  prediction  (Fort  Irwin,  CA) 

To  illustrate  the  issues  associated  with  simulations  for  large  areas  we  use  an  example  of  a  mountain¬ 
ous  region  at  Fort  Irwin,  California.  The  standard  30m  DEM  available  for  the  entire  study  area  (3000  square 
km)  represents  4  million  grid  cells,  a  challenging  data  set  for  currently  available  process-based  simulation 
tools  and  workstations,  at  the  resolution  hardly  sufficient  for  rough  identification  of  high  erosion  risk  areas. 
The  DEM  at  5m  resolution  required  for  more  detailed  modeling  efforts  would  produce  121  million  grid  cells. 
Simulations  would  be  prohibitively  expensive,  if  not  impossible,  with  the  current  computational  resources.  It 
is  clear,  that  for  such  a  large  area,  modeling  at  different  scales  and  resolutions  is  needed,  depending  on  the 
importance  and  complexity  of  the  watersheds.  It  is  important  to  note  that  our  aim  in  this  application  is  to 


Figure  14.  Impact  of  land  use  change  on  erosion  and  deposition:  (a)  original  land  use,  (b)  net  erosion/deposition  for  the 
original  land  use,  (c)  proposed  land  use  change  including  clearing,  thinning,  revegetation  and  stream  buffers,  and  (d) 
USPED  predicted  net  erosion/deposition  for  the  proposed  land  use  showing  20x  increase  in  net  erosion.  Without  the 
protective  measures,  the  increase  in  net  erosion  would  be  33x. 
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illustrate  the  possibility  of  using  standard  elevation  data  for  erosion  simulations  in  a  large  area.  We  do  not 
take  into  account  the  variability  in  climate,  soil  and  cover  properties. 

To  demonstrate  the  impact  of  resolution,  noise  and  systematic  errors  in  standard  elevation  data,  we 
computed  tangential  curvature  for  a  part  of  Fort  Irwin  using  the  standard  30m  DEM.  The  tangential  curvature 
map  was  then  draped  over  a  terrain  map.  The  tangential  curvature  showed  acceptable  structure  in  steeper 
mountainous  areas  but  significant  noise  and  systematic  errors  (stripes)  in  lowland  areas  (Figure  1 5a).  After 
smoothing  and  resampling  to  10m  resolution  using  the  RST  interpolation  method,  the  noise  was  reduced  and 
the  major  topographic  features  became  more  visible  (Figure  15b).  The  analysis  clearly  demonstrates  that  the 
need  for  precision  and  accuracy  in  elevation  data  is  spatially  variable.  Areas  of  flat  terrain  are  much  more 
sensitive  to  noise  and  systematic  errors  in  elevation  data  than  mountains. 

Topographic  potential  for  detachment  limited  erosion  and  transport  capacity  limited  erosion/deposi¬ 
tion  were  estimated  by  the  modified  USLE  (see  Section  4.1)  and  the  USPED  models  (see  Section  4.2), 
respectively.  To  illustrate  the  impact  of  smoothing  and  resampling  on  erosion  modeling,  we  computed  the 
erosion  estimates  using  DEMs  at  90,  30  and  10m.  The  90  and  1 0  m  resolution  DEMs  were  created  by 
resampling  the  original  30  m  DEM.  To  demonstrate  the  differences  in  detail  that  can  be  achieved  at  various 
resolutions,  we  display  the  results  as  color  maps  draped  over  the  10m  resolution  DEM  for  a  36  square  km 
area  (Figures  16  and  17).  The  results  show  that  the  90m  resolution  DEM  is  inadequate  for  capturing  the 
erosion  by  concentrated  flow  with  the  enhanced  USLE  (Figure  16)  and  does  not  allow  the  prediction  of 
realistic  net  erosion/deposition  patterns  using  USPED  (Figure  17).  The  30m  resolution  DEM  starts  to  reveal 
the  structure  of  erosion  patterns,  including  the  concentrated  flow  erosion  in  the  valleys  (Figure  16).  However, 
it  is  not  adequate  for  predicting  net  erosion/deposition  patterns  (Figure  17).  The  10m  resolution  does  not 
improve  the  accuracy  of  the  original  elevation  model,  but  the  smoothing  reduces  the  noise  in  the  data  and  the 
higher  resolution  allows  better  representation  of  terrain  geometry,  thus  leading  to  more  realistic  results  both 
for  the  detachment  limited  erosion  (Figure  16)  and  net  erosion/deposition  (Figure  17). 

To  test  the  applicability  of  the  SIMWE  model  to  large  areas,  we  computed  the  spatial  distribution  of 
steady  state  water  flow,  sediment  flow  and  net  erosion/deposition  for  a  36  square  km  area  at  Fort  Irwin  using 
the  reinterpolated  and  smoothed  10m  DEM  (Figure  18).  The  water  flow  simulation  captured  the  complex 
pattern  of  overland  water  flow,  including  the  dispersal  flow  features  typical  for  this  area,  which  are  difficult 
to  predict  by  more  traditional  approaches  (Figure  18a).  Sediment  flow  rates  were  estimated  by  solution  of 
continuity  of  mass  equation  predicting  high  sediment  flow  rates  in  centers  of  valleys  with  concentrated  flow 
and  dispersal  of  sediment  flow  with  reduction  of  sediment  loads  in  areas  with  alluvial  fans  (Figure  18b).  Net 
erosion/deposition  rates  were  estimated  using  a  transport  capacity  limiting  regime  due  to  the  prevailing  sandy 
soils  (Figure  18c).  The  simulation  shows  the  formation  of  split  gullies  and  their  disappearance  as  the  terrain 
flattens  and  the  alluvial  fans  are  formed. 

This  application  demonstrates  that  processing  of  original  elevation  data  by  adequate  tools,  such  as  in 
our  case  the  RST  method  and  application  of  a  robust,  process-based  model  can  lead  to  a  more  realistic 
prediction  of  a  complex  erosion/deposition  pattern.  This  application  also  proved  that  the  SIMWE  model  can 
be  applied  to  large  areas,  in  spite  of  its  relatively  complex  formulation  and  computational  demands. 


7.5  Evaluation  of  CASC2D-SED  (Goodwin  Creek  watershed,  MS  and  Ft.  Hood,  TX) 

The  CASC2D-SED  model  was  applied  to  Goodwin  Creek  watershed  approximately  60  miles  north  of 
Memphis,  Tennessee  and  the  and  the  Henson  Creek  watershed  at  Fort  Hood,  Texas.  The  Goodwin  Creek 
watershed  is  extensively  gaged  by  the  Agricultural  Research  Service  (ARS)  as  a  research  watershed  in  the 
areas  of  upland  erosion,  in-stream  sediment  transport,  and  watershed  hydrology.  The  watershed  drains  a  total 
area  of  8.26  square  miles.  Only  13%  of  the  Goodwin  Creek  watershed  is  cultivated;  the  remainder  is  idle 
pasture  or  forestland.  The  predominant  soil  texture  for  Goodwin  Creek  watershed  is  silt  loam  with  a  small 
percent  of  sandy  loam. 
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Figure  15.  Tangential  curvature  draped  over  a  DEM  (a)  derived  from  the  original  30  m  DEM  and  (b)  derived  from  a 
smoothed  and  resampled  DEM. 
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Figure  16.  Topographic  potential  for  detachment  limited  erosion  estimated  by  the  enhanced  USLE  for  an  area  at  Fort 
Irwin,  California  at  resolutions  of  90  m  (a),  30  m  (b)  and  10  m  (c). 
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Figure  17.  Erosion  and  deposition  patterns  estimated  by  the  USPED  model  for  an  area  at  Fort  Irwin,  California  at 
resolutions  of  90  m  (a),  30  m  (b)  and  10  m  (c). 
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Figure  18.  Runoff  depth,  sediment  flow  and  erosion/deposition  rates  predicted  by  the  SIMWE  model  from  a  10  m 
resolution  DEM  smoothed  with  the  regularized  spline  with  tension  (RST)  method. 
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At  the  Goodwin  Creek  watershed,  results  of  the  CASC2D-SED  simulations  were  compared  with 
measured  runoff  and  sediment  discharge  from  the  October  17-18,  1981  storm  event.  The  storm  event  began  at 
9:19  pm  and  had  a  total  rainfall  duration  of  3.5  hours  with  very  little  rainfall  preceding  this  event.  Across  the 
watershed,  the  total  rainfall  for  this  event  varied  from  2.55  to  3.11  inches  with  an  average  value  of  2.85 
inches.  A  comparison  of  the  hydrograph  plots  (Figures  19-21)  show  that  CASC2D-SED  was  able  to  consis¬ 
tently  simulate  the  overall  shape  and  rate  of  rise.  The  time  to  peak  was  simulated  within  3%  at  some  places 
(gage  8  and  gage  5),  but  was  off  by  approximately  1 5%  at  gages  2  and  4.  CASC2D-SED  underestimated  the 
total  volume  of  runoff  by  approximately  20%  across  the  watershed.  A  comparison  of  the  sediment  discharge 
plots  (Figures  22-24)  show  that  CASC2D-SED  was  able  to  predict  upland  erosion  from  the  Goodwin  Creek 
Watershed  within  an  acceptable  range  of  -50%  to  200%  of  the  actual  upland  erosion.  This  range  (-50%  to 
200%)  is  generally  considered  by  sedimentation  engineers  to  be  acceptable  when  comparing  computed 
sediment  yields  versus  actual  sediment  yields. 

The  Henson  Creek  watershed  is  contained  almost  completely  within  the  boundaries  of  Fort  Hood, 
Texas.  Most  of  the  watershed  lies  within  the  artillery  impact  area,  parts  of  which  are  subjected  to  frequent  use 
by  military  vehicles.  Soils  range  from  silty  clay  to  clay  loam;  rock  outcrops  are  common.  Model  application 
to  the  Henson  Creek  Watershed,  Fort  Hood,  Texas,  was  performed  using  a  single  storm  occurring  between 
0130  hrs  on  April  25  to  1800  hrs  on  April  27.  Results  indicate  that  CASC2D-SED  provided  reasonably 
accurate  results  (Figure  25). 

For  both  watersheds,  the  CASC2D-SED  model  provided  remarkably  accurate  estimates  of  runoff 
discharge.  Sediment  discharge  estimates  were  less  accurate,  but  generally  within  acceptable  limits.  The 
primaiy  factor  contributing  to  the  errors  in  the  prediction  of  sediment  discharge  appears  to  be  the  fact  that  the 
CASC2D-SED  is  not  designed  to  predict  bank  failure  and  channel  erosion.  As  much  as  70%  of  the  sediment 
discharged  from  the  Goodwin  Creek  watershed  has  been  attributable  to  those  forms  of  erosion. 


7.6  Design  of  erosion  prevention  measures  (Scheyern  Experimental  Farm,  Ger¬ 
many) 

We  first  applied  SIMWE  to  the  Scheyern  Experimental  Farm  in  Germany  to  evaluate  the  erosional 
consequences  of  the  traditional  land  use  management  scenario  that  was  in  use  when  we  first  began  collabora¬ 
tive  efforts  with  the  farm  prior  to  1993  (Figure  9a).  Approximately  80%  of  the  area  was  tilled,  while  20%  was 
maintained  in  permanent  grass  cover.  This  land  use  resulted  in  severe  erosion  in  the  tilled  areas  when  a  large 
storm  event  occurred  during  the  time  when  the  agricultural  fields  were  bare.  The  field  data  showed  that  the 
areas  with  the  highest  water  and  sediment  flow  potential  were  unprotected  and  that  significant  rilling  and 
gully  formation  occurred.  Simulations  performed  for  these  conditions  confirm  high  potential  for  erosion  on 
convex  parts  of  the  bare  hillslopes  and  high  sediment  flows  in  the  areas  of  concentrated  water  flow. 

In  response  to  the  severe  erosion,  the  farm  implemented  “best  management  practices”  in  1993  to 
reduce  erosion  and  enhance  sustainable  agriculture.  The  practices  consisted  of  planting  wide  grassed  buffers 
at  the  valley  bottom.  The  new  scenario  consisted  of  approximately  40%  grassed  areas,  with  the  remainder 
under  cultivation.  Based  on  SIMWE  estimations,  the  improved  land  design  altered  the  spatial  patterns  of 
erosion  and  deposition  within  the  watershed  and  reduced  the  sediment  load  in  the  stream  (Figure  9b). 

We  next  used  the  SIMWE  model  to  investigate  the  possibility  of  redesigning  the  land  use  patterns  to 
minimize  net  soil  loss  and  sediment  loads.  First,  we  used  the  model  to  identify  locations  with  the  highest 
erosion  risk,  assuming  a  uniform  land  use.  Then,  the  original  20%  protective  grass  cover  was  redistributed  to 
the  highest  risk  areas  (Figure  9c).  We  performed  a  simulation  with  the  new  land  use  scenario  to  evaluate  its 
effectiveness.  The  results  demonstrate  that  the  new  design  has  the  potential  to  dramatically  reduce  soil 
erosion  and  sediment  transport  (Figure  9c).  The  high  sediment  flow  in  the  valley  disappeared  and  was  re- 
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Figure  19.  Actual  and  predicted  runoff  hydrograph  for  October  17-19,  1981  at  gages  1  and  2  in  the  Goodwin  Creek 
watershed,  Mississippi. 
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Figure  20.  Actual  and  predicted  runoff  hydrographs  for  October  17-19,  1981  at  gages  3  and  4  in  the  Goodwin  Creek 
watershed,  Mississippi. 
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Figure  21.  Actual  and  predicted  runoff  hydrographs  for  October  17-19,  1981  at  gages  5  and  8  in  the  Goodwin  Creek 
watershed,  Mississippi. 
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Figure  22.  Actual  and  predicted  sediment  discharge  for  October  17-19,  1981  at  gages  1  and  2  in  the  Goodwin  Creek 
watershed,  Mississippi. 
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Figure  23.  Actual  and  predicted  sediment  discharge  for  October  17-19,  1981  at  gages  3  and  4  in  the  Goodwin  Creek 
watershed,  Mississippi. 
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Figure  24.  Actual  and  predicted  sediment  discharge  for  October  17-18,  1981  at  gages  5  and  8  on  the  Goodwin  Creek 
watershed,  Mississippi. 
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Figure  25.  Actual  and  predicted  runoff  discharge  and  sediment  concentration  for  April  25-27,  1997  for  the  Henson  Creek 
watershed,  Fort  Hood,  Texas. 
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placed  by  light  deposition  within  the  grassed  waterways;  net  erosion  was  significantly  reduced  (Mitas  and 
Mitasova  1998).  It  is  interesting  to  note,  that  the  land  use  design  obtained  by  this  rather  simple  computational 
procedure  was  potentially  at  least  as  effective  in  reducing  erosion  and  sediment  transport  as  a  sustainable  land 
use  design  proposed  and  implemented  at  the  farm  in  1993.  Of  importance  to  the  farm  manager,  the  optimized 
design  removed  less  land  from  cultivation  than  the  “best  management”  scenario  developed  without  the  aid  of 
the  erosion  model. 


7.7  Simulation  of  erosion/deposition  at  the  Combat  Maneuver  Training  Center, 
Germany 

In  order  to  compare  the  results  of  the  enhanced  USLE  and  USPED  models,  we  computed  average 
annual  soil  loss  for  the  entire  installation  using  a  10m  DEM.  The  rainfall-runoff  factor  (R)  was  assumed  to  be 
uniform  across  the  installation.  The  soil  erodibility  factor  (K)  varied  spatially  depending  on  the  soil  type.  The 
cover  factor  (C)  was  based  on  the  1997  landcover  map  at  CMTC.  Given  the  absence  of  agricultural  erosion 
control  practices,  the  conservation  practice  factor  (P)  was  assumed  to  be  1,  such  that  it  had  no  effect  on  the 
resulting  erosion  estimates.  Figures  26  and  27  illustrate  the  distribution  of  soil  erosion  and  deposition  as 
estimated  by  the  USLE  and  USPED  models,  respectively.  The  most  obvious  difference  between  the  maps  is 
the  absence  of  areas  of  sediment  deposition  in  the  USLE  soil  loss  map.  As  noted  previously,  due  to  the  nature 
of  the  equation,  the  USLE  is  incapable  of  predicting  deposition  in  its  standard  format.  The  spatial  distribution 
of  zones  of  high  erosion  (>8 1  ac'1  yr1)  are  similar  for  both  models.  Quantitatively,  the  models  predict  a 
similar  percent  of  the  area  experiencing  higher  levels  of  net  erosion  (Table  1).  However,  the  USPED  model 
results  suggest  that  a  significant  portion  of  the  area  predicted  by  the  USLE  to  experience  light  to  moderate 
erosion  (<8 1  ac'1  yr1)  is  more  likely  experiencing  deposition  of  sediments.  Hence,  the  USLE  will  tend  to 
overpredict  net  erosion. 


Table  1.  Quantitative  comparison  of  the  results  of  erosion  and  deposition  estimation  by  the  USLE  and 
USPED  models  at  Combat  Maneuver  Training  Center,  Hohenfels,  Germany. 


Erosion/Deposition  Estimate 

USLE 

USPED 

>0  t  ac'1  yr'1  deposition 

0% 

25% 

0-4  t  ac'1  yr'1  erosion 

79% 

64% 

4-8  t  ac'1  yr'1  erosion 

12% 

4% 

8-12  t  ac'1  yr"1  erosion 

4% 

2% 

12-16  t  ac'1  yr"1  erosion 

2% 

1% 

>16  t  ac'1  yr'1  erosion 

3% 

5% 

The  SIMWE  model  is  computationally  more  demanding  than  the  USPED  model.  Accordingly,  we 
applied  both  models  at  10m  resolution  for  single  watershed.  For  the  USPED  model,  we  used  the  same  param¬ 
eters  as  discussed  above.  For  the  SIMWE  model,  we  performed  the  simulations  for  the  current  land  cover 
considering  two  different  conditions:  (a)  fully  saturated  soil  with  all  rainfall  contributing  to  surface  runoff, 
and  (b)  spatially  variable  rainfall  excess  with  high  portion  of  rainfall  infiltrated  in  the  undisturbed  areas 
covered  mostly  by  forest  and  with  reduced  infiltration  in  the  disturbed  areas.  The  latter  condition  involves 
more  complex  modeling  with  spatially  variable  rainfall  excess  which  had  not  previously  been  tested.  The 
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Figure  26.  Soil  erosion  estimated  by  the  enhanced  USLE  model  at  the  Combat  Maneuver  Training  Center,  Germany. 


Figure  27.  Soil  erosion/deposition  estimated  by  the  USPED  model  at  the  Combat  Maneuver  Training  Center,  Germany. 
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input  parameters  (elevation  gradient,  rainfall  excess,  erodibility,  transportability,  roughness  [Manning’s  n], 
critical  shear  stress)  were  estimated  using  the  elevation,  C-factor  and  K-factor  data,  as  well  as  the  information 
about  the  land  cover  and  soils  from  the  report  by  Warren  and  Kowalski  1989.  The  modeling  results  for  the 
USPED  (Figure  28)  and  SIMWE  (Figure  29)  models  were  similar,  with  the  notable  exception  that  the 
SIMWE  model  predicted  greater  risk  for  erosion  in  hollows  and  centers  of  valleys  and  on  disturbed  ridges. 
The  SIMWE  model  predicted  severe  erosion  risk  over  5.9%  of  the  watershed  with  uniformly  saturated  soil 
conditions.  Spatially  variable  infiltration  increased  the  difference  between  sediment  loads  from  disturbed  and 
undisturbed  areas,  with  undisturbed  areas  being  much  more  stable  and  producing  very  little  sediment  com¬ 
pared  to  disturbed  areas  (Figure  29).  The  impact  of  higher  infiltration  rates  in  undisturbed  areas  reduced  the 
spatial  extent  of  areas  with  severe  erosion  to  3.9%. 

This  application  reveals  several  issues  relevant  to  land  management:  (a)  the  importance  of  erosion 
prevention  measures  in  areas  of  concentrated  flow  and  on  disturbed  ridges,  (b)  the  reduction  of  actual  sedi¬ 
ment  delivery  to  streams  due  to  deposition,  and  (c)  the  stabilizing  effect  of  high  infiltration  in  undisturbed 
areas. 


Figure  28.  Soil  erosion/deposition  estimated  by  the  USPED  model  for  an  area  in  the  eastern  half  of  the  Combat  Maneu¬ 
ver  Training  Center,  Germany. 
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Figure  29.  Soil  erosion/deposition  estimated  by  the  SIMWE  model  for  an  area  in  the  eastern  half  of  the  Combat  Maneu¬ 
ver  Training  Center,  Germany  assuming  uniformly  saturated  soil  conditions  (a)  and  spatially  variable  infiltration  (b)  based 
on  land  use  patterns. 
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Integration  of  topographic  analysis  and  erosion  modeling  with  a  GIS  provides  an  environment  for 
effective  evaluation  of  various  approaches  to  erosion/deposition  risk  assessment  for  landscape  scale  applica¬ 
tions.  This  SERDP  project  has  produced  three  erosion/deposition  models  that  are  applicable  to  a  wide  range 
of  needs.  The  USPED  model  is  a  significant  enhancement  to  the  USLE.  It  is  relatively  simple  to  use  and 
requires  the  least  amount  of  data.  It  is  a  semi-empirical  model  that  produces  estimates  of  long-term  average 
annual  erosion  and  deposition  in  complex  landscapes.  The  CASC2D-SED  is  an  enhancement  to  the  CASC2D 
rainfall-runoff  model.  It  is  also  semi-empirical,  but  provides  estimates  of  erosion  and  deposition  for  indi¬ 
vidual  rainfall  events.  The  SIMWE  model  is  a  physically-based  model  that  can  provide  estimates  of  erosion, 
sediment  transport  and  deposition  for  single  rainfall  events.  Long-term  estimates  of  erosion  and  deposition 
can  be  generated  by  the  CASC2D-SED  and  SIMWE  models  by  simulating  a  realistic  series  of  rainfall  events. 
The  SIMWE  model  requires  much  more  data  input  than  the  other  models.  In  addition,  it  is  computationally 
demanding.  All  three  models  represent  significant  improvements  to  the  field  of  erosion/deposition  modeling. 
For  any  given  application,  the  choice  of  the  model  be  determined  by  the  proposed  goals  and  the  quantity  and 
quality  of  data  to  drive  the  models. 

During  the  process  of  developing  and  applying  the  USPED,  CASC2D-SED  and  SIMWE  models,  we 
have  made  numerous  contributions  to  the  general  understanding  of  erosion  and  deposition  processes  as  well 
as  the  mechanisms  to  include  these  and  other  phenomena  in  multidimensional  modeling.  These  include: 

►  Curvature  of  the  slope  profile  along  a  downhill  gradient  (profile  curvature)  and  curvature  of  the  slope 
perpendicular  to  the  downhill  gradient  (tangential  curvature)  are  both  important  in  determining  net 
erosion  and  deposition.  Profile  curvature  controls  the  velocity  of  runolf,  while  tangential  curvature 
controls  convergence  and  divergence  of  flow,  thus  affecting  both  the  volume  and  velocity  of  runoff.  Both 
curvatures  must  be  considered  in  areas  of  complex  topography.  We  have  introduced  methods  to  calculate 
both  forms  of  curvature  from  digital  elevation  models. 

►  For  areas  of  complex  topography,  ‘upslope  contributing  area’  is  a  more  appropriate  parameter  than  ‘slope 
length’  for  capturing  the  influence  of  accumulated  flow  of  runoff.  We  have  developed  a  simple,  continu¬ 
ous  equation  for  the  computation  of  ‘upslope  contributing  area’  in  complex  terrain. 

►  The  Universal  Soil  Loss  Equation  predicts  soil  erosion  for  situations  where  erosion  is  controlled  prima¬ 
rily  by  the  ability  of  flowing  water  to  detach  soil  particles  (detachment  limited  erosion).  Field  studies  tend 
to  indicate  that  landscape  patterns  of  erosion/deposition  are  better  related  to  situations  where  erosion  is 
limited  primarily  by  the  ability  of  runoff  to  transport  detached  soil  particles  (transport  limited  erosion). 
Hence,  the  results  of  the  USPED  model  more  closely  represent  real-world  conditions. 

►  Soil  erodibility  or  detachability  (Kd)  is  determined  largely  by  soil  qualities.  It  significantly  influences  the 
spatial  patterns  of  erosion  and  deposition,  but  has  only  minimal  influence  on  the  amount  of  sediment  load 
in  a  stream.  Hence,  the  common  practice  of  using  of  in-stream  sediment  loads  for  making  land  rehabilita¬ 
tion  and  management  decisions  is  questionable. 

►  Soil  transportability  (Kt)  is  determined  by  soil  properties  and  vegetation.  It  influences  both  sediment 
loads  and  the  spatial  patterns  of  erosion  and  deposition. 
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►  Surface  roughness,  while  having  minimal  impact  on  soil  erosion,  significantly  affects  both  the  pattern  and 
magnitude  of  sediment  deposition. 

►  Rainfall  excess,  or  the  rainfall  intensity  minus  infiltration  rate  per  unit  time,  affects  the  magnitude  of 
erosion  and  deposition,  but  has  minimal  impact  on  the  spatial  patterns  of  erosion  and  deposition. 

►  The  ability  of  the  soil  to  resist  the  critical  shearing  forces  of  flowing  water  (critical  shear  stress)  primarily 
affects  the  patterns  of  erosion  and  deposition.  However,  on  steep  slopes  and  in  areas  of  concentrated 
flow,  the  ability  of  the  soil  to  resist  shearing  in  an  upslope  position  may  increase  the  magnitude  of  erosion 
downslope  because  clean  water  has  a  greater  potential  to  suspend  sediment  than  dirty  water. 

►  The  various  parameters  of  the  physically-based  SIMWE  model  do  not  act  independently.  It  is  the  inter¬ 
play  of  the  various  parameters  that  ultimately  determines  the  magnitude  and  spatial  distribution  of  erosion 
and  deposition. 

Based  on  our  applications  and  case  studies,  we  are  able  to  make  the  following  conclusions  and 

recommendations: 

►  For  effective  modeling  of  erosion  and  deposition  in  complex  terrain,  a  digital  elevation  model  with  a 
horizontal  resolution  of  20  m  or  less  and  a  vertical  resolution  of  less  than  1  m  is  required. 

►  For  the  purposes  of  erosion/deposition  modeling,  areas  of  flat  terrain  are  more  susceptible  to  noise  and 
systematic  errors  in  elevation  data  than  steeper  terrain. 

►  New-generation  erosion/deposition  models  can  be  valuable  planning  tools  to  effectively  predict  the 
consequences  of  proposed  land-use  changes. 

►  New-generation  erosion/deposition  models  can  be  used  to  effectively  plan  erosion  and  sediment  control 
practices. 

In  support  of  the  new  generation  of  erosion  and  deposition  models  we  have  developed  a  number  of 

geographic  information  system  tools  for  data  processing  and  graphic  display. 
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Peer  reviewed  publications 

Mitas,  L.  and  H.  Mitasova.  1998.  Distributed  soil  erosion  simulation  for  effective  erosion  prevention.  Water 
Resources  Research  34:505-5 1 6. 
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Ogden,  F.L.,  B.  Saghafian  and  W.  F.  Krajewski.  1994.  GIS-based  channel  extraction  and  smoothing  algorithm 
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Kosinovsky,  I.,  H.  Mitasova  and  D.P.  Gerdes.  1994.  Library  for  Multidimensional  Surface  Modeling  and 
Analysis  in  GRASS.  9th  Annual  GRASS/GIS  Conference,  Reston,  VA. 

Brown,  W.M.  1994.  Multi-Dimensional  Visualization  using  GRASS.  9th  Annual  GRASS/GIS  Conference, 
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Trade  journal  articles: 
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The  following  GRASS  GIS  tools  were  developed  during  the  course  of  this  project.  Please  use  the 
following  web  address  to  access  the  manuals  and  tutorials: 

http://www2.gis.uiuc.edu:2280/modviz/reports/cerl98/aman98.html 

1.  gsurf  library  overview  http://www2.gis.uiuc.edu:2280/modviz/reports/cerl98/gsurf/gsurf.html 

2.  gsurf  library  prototypes  and  selected  include  files 

1.1.  public  function  prototypes  http://www2.gis.uiuc.edu:2280/modviz/reports/cerl98/gsurf/ 
protos.html 

1.2.  public  include  file  gsurf.h  http://www2.gis.uiuc.edu:2280/modviz/reports/cerl98/gsurf/ 
gsurf_h.html 

1.3.  public  include  file  keyframe.h  http://www2.gis.uiuc. edu:2280/modviz/reports/cerl98/gsurf/ 
keyffame_h.html 

1.4.  public  color  packing  utility  macros  rgbpack.h  http://www2.gis.uiuc.edu:2280/modviz/reports/ 
cerl98/gsurf/rgbpack_h.html 

1.5.  private  types  and  defines  gstypes.h  http://www2.gis.uiuc. edu:2280/modviz/reports/cerl98/ 
gsurf/gstypes_h.html 

1 .6.  private  utilities  gsget.h  http://www2.gis.uiuc.edu:2280/modviz/reports/cerl98/gsurf/ 
gsget_h.html 

3.  nviz  tutorial  http://www2.gis.uiuc.edu:2280/modviz/reports/cerl98/nviz_tut.html 

4.  sg3d  updated  features  http://www2.gis.uiuc.edu:2280/modviz/reports/cerl98/sg3d42.html 

5.  xganim  user’s  manual  page  http://www2.gis.uiuc.edu:2280/modviz/reports/cerl98/xganim_man.html 

6.  r.out.mpeg  user’s  manual  page  http://www2.gis.uiuc.edu:2280/modviz/reports/cerl98/ 
rout.mpegman.html 

7.  d.siter  user’s  manual  page  http://www2.gis.uiuc.edu:2280/modviz/reports/cerl98/dsiter.html 

8.  d.sites.qual  user’s  manual  page  http://www2.gis.uiuc.edu:2280/modviz/reports/cerl98/dsites_qual.html 

9.  s.info  user’s  manual  page  http://www2.gis.uiuc.edu:2280/modviz/reports/cerl98/sinfo.html 

10.  p.vrml  user’s  manual  page  http://www2.gis.uiuc.edu:2280/modviz/reports/cerl98/pvrml.html 

11.  r3.mkdspf  user’s  manual  page  http://www2.gis.uiuc.edu:2280/modviz/reports/cerl98/r3 .mkdspf.html 

12.  r3.showdspf  user’s  manual  page  http://www2.gis.uiuc.edu:2280/modviz/reports/cerl98/r3.showdspf.html 

13.  s.surf.rst  user’s  manual  page  http://www2.gis.uiuc.edu:2280/modviz/reports/cerl98/ssurfrst.html 

14.  r.flow  user’s  manual  page  http://www2.gis.uiuc.edu:2280/modviz/reports/cerl98/rflow42.html 

15.  r.slope.aspect  user’s  manual  page  http://www2.gis.uiuc.edu:2280/modviz/reports/cerl98/rslopeasp.html 
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